From 6e500b765c737d272d57cb00df760edf87191b25 Mon Sep 17 00:00:00 2001
From: Torkel <torkel.loman@gmail.com>
Date: Thu, 6 Jun 2024 11:47:34 -0400
Subject: [PATCH 1/4] init

---
 .../serialisation_support.jl                  | 23 ++---
 .../serialise_fields.jl                       | 98 +++++++++----------
 .../serialise_reactionsystem.jl               | 10 +-
 .../reactionsystem_serialisation.jl           | 33 +++++--
 4 files changed, 89 insertions(+), 75 deletions(-)

diff --git a/src/reactionsystem_serialisation/serialisation_support.jl b/src/reactionsystem_serialisation/serialisation_support.jl
index 62b366ca5b..40a90a5f46 100644
--- a/src/reactionsystem_serialisation/serialisation_support.jl
+++ b/src/reactionsystem_serialisation/serialisation_support.jl
@@ -36,7 +36,7 @@ function push_field(file_text::String, rn::ReactionSystem, annotate::Bool, top_l
     has_component, get_comp_string, get_comp_annotation = comp_funcs
     has_component(rn) || (return (file_text, false))
 
-    # Prepares the text creating the field. For non-top level systems, adds `local `. Observables
+    # Prepares the text creating the field. For non-top level systems, add `local `. Observables
     # must be handled differently (as the declaration is not at the beginning of the code for these).
     # The independent variables is not declared as a variable, and also should not have a `1ocal `.
     write_string = get_comp_string(rn)
@@ -54,7 +54,7 @@ function push_field(file_text::String, rn::ReactionSystem, annotate::Bool, top_l
     return (file_text * write_string, true)
 end
 
-# Generic function for creating an string for an unsupported argument.
+# Generic function for creating an string for a unsupported argument.
 function get_unsupported_comp_string(component::String)
     @warn "Writing ReactionSystem models with $(component) is currently not supported. This field is not written to the file."
     return ""
@@ -82,9 +82,10 @@ function syms_2_strings(syms)
     return get_substring_end("$(convert(Vector{Any}, strip_called_syms))", 4, 0)
 end    
 
-# Converts a vector of symbolics (e.g. the species or parameter vectors) to a string corresponding to 
-# the code required to declare them (potential @parameters or @species commands must still be added).
-# The `multiline_format` option formats it with a `begin ... end` block and declarations on separate lines.
+# Converts a vector of symbolic variables (e.g. the species or parameter vectors) to a string
+# corresponding to the code required to declare them (potential @parameters or @species commands
+# must still be added). The `multiline_format` option formats it with a `begin ... end` block
+# and declarations on separate lines.
 function syms_2_declaration_string(syms; multiline_format = false)
     decs_string = (multiline_format ? " begin" : "")
     for sym in syms
@@ -102,7 +103,7 @@ function sym_2_declaration_string(sym; multiline_format = false)
     # Creates the basic symbol. The `"$(sym)"` ensures that we get e.g. "X(t)" and not "X".
     dec_string = "$(sym)"
 
-    # If the symbol have a non-default type, appends the declaration of this.
+    # If the symbol has a non-default type, appends the declaration of this.
     # Assumes that the type is on the form `SymbolicUtils.BasicSymbolic{X}`. Contain error checks
     # to ensure that this is the case.
     if !(sym isa SymbolicUtils.BasicSymbolic{Real})
@@ -136,7 +137,7 @@ end
 
 # Converts a generic value to a String. Handles each type of value separately. Unsupported values might
 # not necessarily generate valid code, and hence throw errors. Primarily used to write default values
-# and metadata values (which hopefully almost exclusively) has simple, supported, types. Ideally,
+# and metadata values (which hopefully almost exclusively) have simple, supported, types. Ideally,
 # more supported types can be added here.
 x_2_string(x::Num) = expression_2_string(x)
 x_2_string(x::SymbolicUtils.BasicSymbolic{<:Real}) = expression_2_string(x)
@@ -261,7 +262,7 @@ function get_dep_syms(sym)
     return Symbolics.get_variables(ModelingToolkit.getdefault(sym))
 end
 
-# Checks if a symbolic depends on an symbolics in a vector being declared.
+# Checks if a symbolic depends on a symbolics in a vector being declared.
 # Because Symbolics has to utilise `isequal`, the `isdisjoint` function cannot be used.
 function depends_on(sym, syms)
     dep_syms = get_dep_syms(sym)
@@ -275,8 +276,8 @@ end
 
 # For a set of remaining parameters/species/variables (remaining_syms), return this split into
 # two sets:
-# One with those that does not depend on any sym in `all_remaining_syms`.
-# One with those that does depend on at least one sym in `all_remaining_syms`.
+# One with those that do not depend on any sym in `all_remaining_syms`.
+# One with those that do depend on at least one sym in `all_remaining_syms`.
 # The first set is returned. Next `remaining_syms` is updated to be the second set.
 function dependency_split!(remaining_syms, all_remaining_syms)
     writable_syms = filter(sym -> !depends_on(sym, all_remaining_syms), remaining_syms)
@@ -289,7 +290,7 @@ end
 
 
 # Checks if a symbolic's declaration is "complicated". The declaration is considered complicated
-# if it have metadata, default value, or type designation that must be declared.
+# if it has metadata, default value, or type designation that must be declared.
 function complicated_declaration(sym)
     isempty(get_metadata_to_declare(sym)) || (return true)
     ModelingToolkit.hasdefault(sym) && (return true)
diff --git a/src/reactionsystem_serialisation/serialise_fields.jl b/src/reactionsystem_serialisation/serialise_fields.jl
index 964e3bb617..09ae2fe66a 100644
--- a/src/reactionsystem_serialisation/serialise_fields.jl
+++ b/src/reactionsystem_serialisation/serialise_fields.jl
@@ -1,7 +1,7 @@
 ### Handles Independent Variables ###
 
-# Checks if the reaction system have any independent variable. True for all valid reaction systems.
-function has_iv(rn::ReactionSystem)
+# Checks if the reaction system has any independent variable. True for all valid reaction systems.
+function seri_has_iv(rn::ReactionSystem)
     return true
 end
 
@@ -17,13 +17,13 @@ function get_iv_annotation(rn::ReactionSystem)
 end
 
 # Combines the 3 independent variable-related functions in a constant tuple.
-IV_FS = (has_iv, get_iv_string, get_iv_annotation)
+IV_FS = (seri_has_iv, get_iv_string, get_iv_annotation)
 
 
 ### Handles Spatial Independent Variables ###
 
-# Checks if the reaction system have any spatial independent variables.
-function has_sivs(rn::ReactionSystem)
+# Checks if the reaction system has any spatial independent variables.
+function seri_has_sivs(rn::ReactionSystem)
     return !isempty(get_sivs(rn))
 end
 
@@ -38,7 +38,7 @@ function get_sivs_annotation(rn::ReactionSystem)
 end
 
 # Combines the 3 independent variables-related functions in a constant tuple.
-SIVS_FS = (has_sivs, get_sivs_string, get_sivs_annotation)
+SIVS_FS = (seri_has_sivs, get_sivs_string, get_sivs_annotation)
 
 
 ### Handles Species, Variables, and Parameters ###
@@ -46,15 +46,15 @@ SIVS_FS = (has_sivs, get_sivs_string, get_sivs_annotation)
 # Function which handles the addition of species, variable, and parameter declarations to the file
 # text. These must be handled as a unity in case there are default value dependencies between these.
 function handle_us_n_ps(file_text::String, rn::ReactionSystem, annotate::Bool, top_level::Bool)
-    # Fetches the systems parameters, species, and variables. Computes the `has_` `Bool`s.
+    # Fetches the system's parameters, species, and variables. Computes the `has_` `Bool`s.
     ps_all = get_ps(rn)
     sps_all = get_species(rn)
     vars_all = filter(!isspecies, get_unknowns(rn))
-    has_ps = has_parameters(rn)
-    has_sps = has_species(rn)
-    has_vars = has_variables(rn)
+    has_ps = seri_has_parameters(rn)
+    has_sps = seri_has_species(rn)
+    has_vars = seri_has_variables(rn)
 
-    # Checks which sets have dependencies which requires managing.
+    # Checks which sets have dependencies which require managing.
     p_deps = any(depends_on(p, [ps_all; sps_all; vars_all]) for p in ps_all)
     sp_deps = any(depends_on(sp, [sps_all; vars_all]) for sp in sps_all)
     var_deps = any(depends_on(var, vars_all) for var in vars_all)
@@ -93,7 +93,7 @@ function handle_us_n_ps(file_text::String, rn::ReactionSystem, annotate::Bool, t
 
         # Pre-declares the sets with written/remaining parameters/species/variables.
         # Whenever all/none are written depends on whether there were any initial dependencies.
-        # `deepcopy` is required as these gets mutated by `dependency_split!`. 
+        # `deepcopy` is required as these get mutated by `dependency_split!`. 
         remaining_ps = (p_deps ? deepcopy(ps_all) : [])
         remaining_sps = (sp_deps ? deepcopy(sps_all) : [])
         remaining_vars = (var_deps ? deepcopy(vars_all) : [])
@@ -113,7 +113,7 @@ function handle_us_n_ps(file_text::String, rn::ReactionSystem, annotate::Bool, t
             isempty(writable_vars) || @string_append! us_n_ps_string get_variables_string(writable_vars) "\n"
         end
 
-        # For parameters, species, and/or variables with dependencies, creates final vectors.
+        # For parameters, species, and/or variables with dependencies, create final vectors.
         p_deps && (@string_append! us_n_ps_string "ps = " syms_2_strings(ps_all) "\n")
         sp_deps && (@string_append! us_n_ps_string "sps = " syms_2_strings(sps_all) "\n")
         var_deps && (@string_append! us_n_ps_string "vars = " syms_2_strings(vars_all) "\n")
@@ -127,17 +127,17 @@ function handle_us_n_ps(file_text::String, rn::ReactionSystem, annotate::Bool, t
         us_n_ps_string = replace(us_n_ps_string, "\nvars = " => "\nlocal vars = ")
     end
 
-    # Merges the file text with `us_n_ps_string` and return the final outputs.
+    # Merges the file text with `us_n_ps_string` and returns the final outputs.
     return file_text * us_n_ps_string, has_ps, has_sps, has_vars
 end
 
 
 ### Handles Parameters ###
-# Unlike most other fields, there are not called via `push_field`, but rather via `handle_us_n_ps`.
+# Unlike most other fields, these are not called via `push_field`, but rather via `handle_us_n_ps`.
 # Hence they work slightly differently.
 
-# Checks if the reaction system have any parameters.
-function has_parameters(rn::ReactionSystem)
+# Checks if the reaction system has any parameters.
+function seri_has_parameters(rn::ReactionSystem)
     return !isempty(get_ps(rn))
 end
 
@@ -156,11 +156,11 @@ end
 
 
 ### Handles Species ###
-# Unlike most other fields, there are not called via `push_field`, but rather via `handle_us_n_ps`.
+# Unlike most other fields, these are not called via `push_field`, but rather via `handle_us_n_ps`.
 # Hence they work slightly differently.
 
-# Checks if the reaction system have any species.
-function has_species(rn::ReactionSystem)
+# Checks if the reaction system has any species.
+function seri_has_species(rn::ReactionSystem)
     return !isempty(get_species(rn))
 end
 
@@ -179,11 +179,11 @@ end
 
 
 ### Handles Variables ###
-# Unlike most other fields, there are not called via `push_field`, but rather via `handle_us_n_ps`.
+# Unlike most other fields, these are not called via `push_field`, but rather via `handle_us_n_ps`.
 # Hence they work slightly differently.
 
-# Checks if the reaction system have any variables.
-function has_variables(rn::ReactionSystem)
+# Checks if the reaction system has any variables.
+function seri_has_variables(rn::ReactionSystem)
     return length(get_unknowns(rn)) > length(get_species(rn))
 end
 
@@ -201,13 +201,13 @@ function get_variables_annotation(rn::ReactionSystem)
 end
 
 # Combines the 3 variables-related functions in a constant tuple.
-VARIABLES_FS = (has_variables, get_variables_string, get_variables_annotation)
+VARIABLES_FS = (seri_has_variables, get_variables_string, get_variables_annotation)
 
 
 ### Handles Reactions ###
 
-# Checks if the reaction system have any reactions.
-function has_reactions(rn::ReactionSystem)
+# Checks if the reaction system has any reactions.
+function seri_has_reactions(rn::ReactionSystem)
     return length(reactions(rn)) != 0
 end
 
@@ -265,14 +265,14 @@ function get_reactions_annotation(rn::ReactionSystem)
     return "Reactions:"
 end
 
-# Combines the 3 reactions-related functions in a constant tuple.
-REACTIONS_FS = (has_reactions, get_reactions_string, get_reactions_annotation)
+# Combines the 3 reaction-related functions in a constant tuple.
+REACTIONS_FS = (seri_has_reactions, get_reactions_string, get_reactions_annotation)
 
 
 ### Handles Equations ###
 
-# Checks if the reaction system have any equations.
-function has_equations(rn::ReactionSystem)
+# Checks if the reaction system has any equations.
+function seri_has_equations(rn::ReactionSystem)
     return length(get_eqs(rn)) > length(get_rxs(rn))
 end
 
@@ -302,13 +302,13 @@ function get_equations_annotation(rn::ReactionSystem)
 end
 
 # Combines the 3 equations-related functions in a constant tuple.
-EQUATIONS_FS = (has_equations, get_equations_string, get_equations_annotation)
+EQUATIONS_FS = (seri_has_equations, get_equations_string, get_equations_annotation)
 
 
 ### Handles Observables ###
 
-# Checks if the reaction system have any observables.
-function has_observed(rn::ReactionSystem)
+# Checks if the reaction system has any observables.
+function seri_has_observed(rn::ReactionSystem)
     return !isempty(observed(rn))
 end
 
@@ -353,13 +353,13 @@ function get_observed_annotation(rn::ReactionSystem)
 end
 
 # Combines the 3 -related functions in a constant tuple.
-OBSERVED_FS = (has_observed, get_observed_string, get_observed_annotation)
+OBSERVED_FS = (seri_has_observed, get_observed_string, get_observed_annotation)
 
 
 ### Handles Continuous Events ###
 
-# Checks if the reaction system have any continuous events.
-function has_continuous_events(rn::ReactionSystem)
+# Checks if the reaction system have has continuous events.
+function seri_has_continuous_events(rn::ReactionSystem)
     return !isempty(MT.get_continuous_events(rn))
 end
 
@@ -410,13 +410,13 @@ function get_continuous_events_annotation(rn::ReactionSystem)
 end
 
 # Combines the 3 -related functions in a constant tuple.
-CONTINUOUS_EVENTS_FS = (has_continuous_events, get_continuous_events_string, get_continuous_events_annotation)
+CONTINUOUS_EVENTS_FS = (seri_has_continuous_events, get_continuous_events_string, get_continuous_events_annotation)
 
 
 ### Handles Discrete Events ###
 
-# Checks if the reaction system have any discrete events.
-function has_discrete_events(rn::ReactionSystem)
+# Checks if the reaction system has any discrete events.
+function seri_has_discrete_events(rn::ReactionSystem)
     return !isempty(MT.get_discrete_events(rn))
 end
 
@@ -466,17 +466,17 @@ function get_discrete_events_annotation(rn::ReactionSystem)
 end
 
 # Combines the 3 -related functions in a constant tuple.
-DISCRETE_EVENTS_FS = (has_discrete_events, get_discrete_events_string, get_discrete_events_annotation)
+DISCRETE_EVENTS_FS = (seri_has_discrete_events, get_discrete_events_string, get_discrete_events_annotation)
 
 
 ### Handles Systems ###
 
 # Specific `push_field` function, which is used for the system field (where the annotation option
 # must be passed to the `get_component_string` function). Since non-ReactionSystem systems cannot be 
-# written to file, this functions throws an error if any such systems are encountered.
+# written to file, this function throws an error if any such systems are encountered.
 function push_systems_field(file_text::String, rn::ReactionSystem, annotate::Bool, top_level::Bool)
-    # Checks whther there are any subsystems, and if these are ReactionSystems.
-    has_systems(rn) || (return (file_text, false))
+    # Checks whether there are any subsystems, and if these are ReactionSystems.
+    seri_has_systems(rn) || (return (file_text, false))
     if any(!(system isa ReactionSystem) for system in MT.get_systems(rn)) 
         error("Tries to write a ReactionSystem to file which have non-ReactionSystem subs-systems. This is currently not possible.")
     end
@@ -489,8 +489,8 @@ function push_systems_field(file_text::String, rn::ReactionSystem, annotate::Boo
     return (file_text * write_string, true)
 end
 
-# Checks if the reaction system have any systems.
-function has_systems(rn::ReactionSystem)
+# Checks if the reaction system has any systems.
+function seri_has_systems(rn::ReactionSystem)
     return !isempty(MT.get_systems(rn))
 end
 
@@ -519,13 +519,13 @@ function get_systems_annotation(rn::ReactionSystem)
 end
 
 # Combines the 3 systems-related functions in a constant tuple.
-SYSTEMS_FS = (has_systems, get_systems_string, get_systems_annotation)
+SYSTEMS_FS = (seri_has_systems, get_systems_string, get_systems_annotation)
 
 
 ### Handles Connection Types ###
 
-# Checks if the reaction system have any connection types.
-function has_connection_type(rn::ReactionSystem)
+# Checks if the reaction system has any connection types.
+function seri_has_connection_type(rn::ReactionSystem)
     return false
 end
 
@@ -540,4 +540,4 @@ function get_connection_type_annotation(rn::ReactionSystem)
 end
 
 # Combines the 3 connection types-related functions in a constant tuple.
-CONNECTION_TYPE_FS = (has_connection_type, get_connection_type_string, get_connection_type_annotation)
\ No newline at end of file
+CONNECTION_TYPE_FS = (seri_has_connection_type, get_connection_type_string, get_connection_type_annotation)
\ No newline at end of file
diff --git a/src/reactionsystem_serialisation/serialise_reactionsystem.jl b/src/reactionsystem_serialisation/serialise_reactionsystem.jl
index 3dc61c5c75..7885c6a105 100644
--- a/src/reactionsystem_serialisation/serialise_reactionsystem.jl
+++ b/src/reactionsystem_serialisation/serialise_reactionsystem.jl
@@ -53,9 +53,9 @@ function get_full_system_string(rn::ReactionSystem, annotate::Bool, top_level::B
     file_text = ""
     
     # Goes through each type of system component, potentially adding it to the string.
-    # Species, variables, and parameters must be handled differently in case there is default-values
+    # Species, variables, and parameters must be handled differently in case there are default values
     # dependencies between them. 
-    # Systems uses custom `push_field` function as these requires the annotation `Bool`to be passed 
+    # Systems use custom `push_field` function as these require the annotation `Bool`to be passed 
     # to the function that creates the next sub-system declarations.
     file_text, _ = push_field(file_text, rn, annotate, top_level, IV_FS)
     file_text, has_sivs = push_field(file_text, rn, annotate, top_level, SIVS_FS)
@@ -68,8 +68,8 @@ function get_full_system_string(rn::ReactionSystem, annotate::Bool, top_level::B
     file_text, has_systems = push_systems_field(file_text, rn, annotate, top_level)
     file_text, has_connection_type = push_field(file_text, rn, annotate, top_level, CONNECTION_TYPE_FS)
 
-    # Finalises the system. Creates the final `ReactionSystem` call.
-    # Enclose everything ing a `let ... end` block.
+    # Finalise the system. Creates the final `ReactionSystem` call.
+    # Enclose everything in a `let ... end` block.
     rs_creation_code = make_reaction_system_call(rn, annotate, top_level, has_sivs, has_species, 
                                                  has_variables, has_parameters, has_reactions, 
                                                  has_equations, has_observed, has_continuous_events, 
@@ -82,7 +82,7 @@ function get_full_system_string(rn::ReactionSystem, annotate::Bool, top_level::B
 end
 
 # Creates a ReactionSystem call for creating the model. Adds all the correct inputs to it. The input
-# `has_` `Bool`s described which inputs are used. If model is `complete`, this is handled here.
+# `has_` `Bool`s described which inputs are used. If the model is `complete`, this is handled here.
 function make_reaction_system_call(rs::ReactionSystem, annotate, top_level, has_sivs, has_species, 
                                    has_variables, has_parameters, has_reactions, has_equations, 
                                    has_observed, has_continuous_events, has_discrete_events,
diff --git a/test/miscellaneous_tests/reactionsystem_serialisation.jl b/test/miscellaneous_tests/reactionsystem_serialisation.jl
index eab88d0abf..b06c434892 100644
--- a/test/miscellaneous_tests/reactionsystem_serialisation.jl
+++ b/test/miscellaneous_tests/reactionsystem_serialisation.jl
@@ -359,17 +359,11 @@ end
 # Tests for (slightly more) complicate system created via the DSL.
 # Tests for cases where the number of input is untested (i.e. multiple observables and continuous
 # events, but single equations and discrete events).
-# Currently broken due to Symbolics doing something weird with observable variables, where these
-# end up not being internal due to something internal in symbolics. I have tried tracking down the
-# obscure symbolics subfields.
+# Tests with and without `safety_check`.
 let 
     # Declares the model.
     rs = @reaction_network begin
         @equations D(V) ~ 1 - V
-        @observables begin
-            X2 ~ 2*X
-            X3 ~ 3*X
-        end
         @continuous_events begin
             [X ~ 5.0] => [X ~ X + 1.0]
             [X ~ 20.0] => [X ~ X - 1.0]
@@ -379,11 +373,30 @@ let
     end
 
     # Checks that serialisation works.
-    save_reactionsystem("serialised_rs.jl", rs; safety_check = false)
-    @test_broken isequal(rs, include("../serialised_rs.jl"))
-    rm("serialised_rs.jl")
+    save_reactionsystem("serialised_rs_1.jl", rs)
+    save_reactionsystem("serialised_rs_2.jl", rs; safety_check = false)
+    isequal(rs, include("../serialised_rs_1.jl"))
+    isequal(rs, include("../serialised_rs_2.jl"))
+    rm("serialised_rs_1.jl")
+    rm("serialised_rs_2.jl")
 end
 
+# Tests for system where species depends on multiple independent variables.
+# Tests for system where variables depends on multiple independent variables.
+let
+    rs = @reaction_network begin
+        @ivs t x y z
+        @parameters p
+        @species X(t,x,y) Y(t,x,y) XY(t,x,y) Z(t,x,y)
+        @variables V(t,x,z)
+        (kB,kD), X + Y <--> XY
+    end
+    save_reactionsystem("serialised_rs.jl", rs)
+    @test !ModelingToolkit.isequal(rs, include("../serialised_rs.jl"))
+    rm("serialised_rs.jl.jl")
+end
+
+
 ### Other Tests ###
 
 # Tests that an error is generated when non-`ReactionSystem` subs-systems are used.

From bb2c293f701cf5c428f6bf02dcdbe3e40337ff33 Mon Sep 17 00:00:00 2001
From: Torkel <torkel.loman@gmail.com>
Date: Thu, 6 Jun 2024 16:18:12 -0400
Subject: [PATCH 2/4] up

---
 test/miscellaneous_tests/reactionsystem_serialisation.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/test/miscellaneous_tests/reactionsystem_serialisation.jl b/test/miscellaneous_tests/reactionsystem_serialisation.jl
index b06c434892..1eac2e6985 100644
--- a/test/miscellaneous_tests/reactionsystem_serialisation.jl
+++ b/test/miscellaneous_tests/reactionsystem_serialisation.jl
@@ -392,7 +392,7 @@ let
         (kB,kD), X + Y <--> XY
     end
     save_reactionsystem("serialised_rs.jl", rs)
-    @test !ModelingToolkit.isequal(rs, include("../serialised_rs.jl"))
+    @test ModelingToolkit.isequal(rs, include("../serialised_rs.jl"))
     rm("serialised_rs.jl.jl")
 end
 

From d3b80d6e7f8cde70f911e365e41ff59ee0eb179a Mon Sep 17 00:00:00 2001
From: Torkel <torkel.loman@gmail.com>
Date: Fri, 7 Jun 2024 10:04:01 -0400
Subject: [PATCH 3/4] writing fix

---
 test/miscellaneous_tests/reactionsystem_serialisation.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/test/miscellaneous_tests/reactionsystem_serialisation.jl b/test/miscellaneous_tests/reactionsystem_serialisation.jl
index 1eac2e6985..bc3845573a 100644
--- a/test/miscellaneous_tests/reactionsystem_serialisation.jl
+++ b/test/miscellaneous_tests/reactionsystem_serialisation.jl
@@ -393,7 +393,7 @@ let
     end
     save_reactionsystem("serialised_rs.jl", rs)
     @test ModelingToolkit.isequal(rs, include("../serialised_rs.jl"))
-    rm("serialised_rs.jl.jl")
+    rm("serialised_rs.jl")
 end
 
 

From 99a9753c32ed887e5237534d6d612be2c71ef187 Mon Sep 17 00:00:00 2001
From: Torkel <torkel.loman@gmail.com>
Date: Fri, 7 Jun 2024 12:12:51 -0400
Subject: [PATCH 4/4] init

---
 src/reaction.jl                      |  6 ++
 test/reactionsystem_core/reaction.jl | 87 ++++++++++++++++++++++++++--
 2 files changed, 89 insertions(+), 4 deletions(-)

diff --git a/src/reaction.jl b/src/reaction.jl
index 23538214aa..7f82f1376d 100644
--- a/src/reaction.jl
+++ b/src/reaction.jl
@@ -163,6 +163,12 @@ end
 function Reaction(rate, subs, prods, substoich, prodstoich;
                   netstoich = nothing, metadata = Pair{Symbol, Any}[], 
                   only_use_rate = metadata_only_use_rate_check(metadata), kwargs...)
+    # Handles empty/nothing vectors.
+    isnothing(subs) || isempty(subs) && (subs = nothing)
+    isnothing(prods) || isempty(prods) && (prods = nothing)
+    isnothing(substoich) || isempty(substoich) && (substoich = nothing)
+    isnothing(prodstoich) || isempty(prodstoich) && (prodstoich = nothing)
+
     (isnothing(prods) && isnothing(subs)) &&
         throw(ArgumentError("A reaction requires a non-nothing substrate or product vector."))
     (isnothing(prodstoich) && isnothing(substoich)) &&
diff --git a/test/reactionsystem_core/reaction.jl b/test/reactionsystem_core/reaction.jl
index 4760641c4b..e70544d0fb 100644
--- a/test/reactionsystem_core/reaction.jl
+++ b/test/reactionsystem_core/reaction.jl
@@ -8,6 +8,89 @@ using ModelingToolkit: value, get_variables!
 # Sets the default `t` to use.
 t = default_t()
 
+### Reaction Constructor Tests ###
+
+# Checks that `Reaction`s can be successfully created using various complicated inputs.
+# Checks that the `Reaction`s have the correct type, and the correct net stoichiometries are generated.
+let
+    # Declare symbolic variables.
+    @parameters k n1 n2::Int32 x [isconstantspecies=true]
+    @species X(t) Y(t) Z(t)
+    @variables A(t)
+
+    # Tries for different types of rates (should not matter).
+    for rate in (k, k*A, 2, 3.0, 4//3)
+        # Creates `Reaction`s.
+        rx1 = Reaction(rate, [X], [])
+        rx2 = Reaction(rate, [x], [Y], [1.5], [1])
+        rx3 = Reaction(rate, [x, X], [], [n1 + n2, n2], [])
+        rx4 = Reaction(rate, [X, Y], [X, Y, Z], [2//3, 3], [1//3, 1, 2])
+        rx5 = Reaction(rate, [X, Y], [X, Y, Z], [2, 3], [1, n1, n2])
+        rx6 = Reaction(rate, [X], [x], [n1], [1])
+
+        # Check `Reaction` types.
+        @test rx1 isa Reaction{Any,Int64}
+        @test rx2 isa Reaction{Any,Float64}
+        @test rx3 isa Reaction{Any,Any}
+        @test rx4 isa Reaction{Any,Rational{Int64}}
+        @test rx5 isa Reaction{Any,Any}
+        @test rx6 isa Reaction{Any,Any}
+
+        # Check `Reaction` net stoichiometries.
+        issetequal(rx1.netstoich, [X => -1])
+        issetequal(rx2.netstoich, [x => -1.5, Y => 1.0])
+        issetequal(rx3.netstoich, [x => -n1 - n2, X => -n2])
+        issetequal(rx4.netstoich, [X => -1//3, Y => -2//1, Z => 2//1])
+        issetequal(rx5.netstoich, [X => -1, Y => n1 - 3, Z => n2])
+        issetequal(rx6.netstoich, [X => -n1, x => 1])
+    end
+end
+
+# Tests that various `Reaction` constructors gives identical inputs.
+let
+    # Declare symbolic variables.
+    @parameters k n1 n2::Int32
+    @species X(t) Y(t) Z(t)
+    @variables A(t)
+
+    # Tests that the three-argument constructor generates correct result.
+    @test Reaction(k*A, [X], [Y, Z]) == Reaction(k*A, [X], [Y, Z], [1], [1, 1])
+
+    # Tests that `[]` and `nothing` can be used interchangeably.
+    @test Reaction(k*A, [X, Z], nothing) == Reaction(k*A, [X, Z], [])
+    @test Reaction(k*A, nothing, [Y, Z]) == Reaction(k*A, [], [Y, Z])
+    @test Reaction(k*A, [X, Z], nothing, [n1 + n2, 2], nothing) == Reaction(k*A, [X, Z], [], [n1 + n2, 2], [])
+    @test Reaction(k*A, nothing, [Y, Z], nothing, [n1 + n2, 2]) == Reaction(k*A, [], [Y, Z], [], [n1 + n2, 2])
+end
+
+# Tests that various incorrect inputs yields errors.
+let
+    # Declare symbolic variables.
+    @parameters k n1 n2::Int32
+    @species X(t) Y(t) Z(t)
+    @variables A(t)
+    
+    # Neither substrates nor products.
+    @test_throws ArgumentError Reaction(k*A, [], [])
+
+    # Substrate vector not of equal length to substrate stoichiometry vector.
+    @test_throws ArgumentError Reaction(k*A, [X, X, Z], [], [1, 2], [])
+
+    # Product vector not of equal length to product stoichiometry vector.
+    @test_throws ArgumentError Reaction(k*A, [], [X, X, Z], [], [1, 2])
+
+    # Repeated substrates.
+    @test_throws ArgumentError Reaction(k*A, [X, X, Z], [])
+
+    # Repeated products.
+    @test_throws ArgumentError Reaction(k*A, [], [Y, Z, Z])
+
+    # Non-valid reactants (parameter or variable).
+    @test_throws ArgumentError Reaction(k*A, [], [A])
+    @test_throws ArgumentError Reaction(k*A, [], [k])
+end
+
+
 ### Test Basic Accessors ###
 
 # Tests the `get_variables` function.
@@ -42,7 +125,6 @@ end
 # Tests basic accessor functions.
 # Tests that repeated metadata entries are not permitted.
 let
-    @variables t
     @parameters k
     @species X(t) X2(t)
 
@@ -60,7 +142,6 @@ end
 
 # Tests accessors for system without metadata.
 let
-    @variables t
     @parameters k
     @species X(t) X2(t)
 
@@ -77,7 +158,6 @@ end
 # Tests basic accessor functions.
 # Tests various metadata types.
 let
-    @variables t
     @parameters k
     @species X(t) X2(t)
 
@@ -109,7 +189,6 @@ end
 
 # Tests the noise scaling metadata.
 let
-    @variables t
     @parameters k η  
     @species X(t) X2(t)