Serialise ReactionSystem to file.

[TorkelE]

Adds a function, save_reaction_network that saves a reaction network to a file:

    save_reaction_network(filename::String, rn::ReactionSystem; annotate = true, safety_check = true)

Save a `ReactionSystem` model to a file. The `ReactionSystem` is saved as runnable Julia code. This
can both be used to save a `ReactionSystem` model, but also to write it to a file for easy inspection.

Arguments:
- `filename`: The name of the file to which the `ReactionSystem` is saved.
- `rn`: The `ReactionSystem` which should be saved to a file.
- `annotate = true`: Whether annotation should be added to the file.
- `safety_check = true`: After serialisation, Catalyst will automatically load the serialised
  `ReactionSystem` and check that it is equal to `rn`. If it is not, an error will be thrown. For
  models without the `connection_type` field, this should not happen. If performance is required 
  (i.e. when saving a large number of models), this can be disabled by setting `safety_check = false`.

It takes all fields into account, except for the connection field (in which case there is a warning).

As an example, the following model:

using Catalyst

# Sets the default `t` and `D` to use.
t = default_t()
D = default_time_deriv()

# Prepares spatial independent variables (technically not used and only supplied to systems).
sivs = @variables x y z [description="A spatial independent variable."]

# Prepares parameters, species, and variables.
@parameters p d k1_1 k2_1 k1_2 k2_2 k1_3 k2_3 k1_4 k2_4 a b_1 b_2 b_3 b_4
@parameters begin
    t_1 = 2.0
    t_2::Float64
    t_3, [description="A parameter."]
    t_4::Float32 = p, [description="A parameter."]
end 
@species X(t) X2_1(t) X2_2(t) X2_3(t) X2_4(t)=p [description="A species."]
@variables A(t)=p [description="A variable."] B_1(t) B_2(t) B_3(t) B_4(t) O_1(t) O_2(t) O_3(t) O_4(t)

# Prepares all equations.
eqs_1 = [
    Reaction(p, [], [X]; metadata = [:description => "A reaction"]),
    Reaction(d, [X], []),
    Reaction(k1_1, [X], [X2_1], [2], [1]),
    Reaction(k2_1, [X2_1], [X], [1], [2]),
    D(A) ~ a - A,
    A + 2B_1^3 ~ b_1 * X
]
eqs_2 = [
    Reaction(p, [], [X]; metadata = [:description => "A reaction"]),
    Reaction(d, [X], []),
    Reaction(k1_2, [X], [X2_2], [2], [1]),
    Reaction(k2_2, [X2_2], [X], [1], [2]),
    D(A) ~ a - A,
    A + 2B_2^3 ~ b_2 * X
]
eqs_3 = [
    Reaction(p, [], [X]; metadata = [:description => "A reaction"]),
    Reaction(d, [X], []),
    Reaction(k1_3, [X], [X2_3], [2], [1]),
    Reaction(k2_3, [X2_3], [X], [1], [2]),
    D(A) ~ a - A,
    A + 2B_3^3 ~ b_3 * X
]
eqs_4 = [
    Reaction(p, [], [X]; metadata = [:description => "A reaction"]),
    Reaction(d, [X], []),
    Reaction(k1_4, [X], [X2_4], [2], [1]),
    Reaction(k2_4, [X2_4], [X], [1], [2]),
    D(A) ~ a - A,
    A + 2B_4^3 ~ b_4 * X
]

# Prepares all observables. 
observed_1 = [O_1 ~ X + 2*X2_1]
observed_2 = [O_2 ~ X + 2*X2_2]
observed_3 = [O_3 ~ X + 2*X2_3]
observed_4 = [O_4 ~ X + 2*X2_4]

# Prepares all events.
continuous_events_1 = [(A ~ t_1) => [A ~ A + 2.0, X ~ X/2]]
continuous_events_2 = [(A ~ t_2) => [A ~ A + 2.0, X ~ X/2]]
continuous_events_3 = [(A ~ t_3) => [A ~ A + 2.0, X ~ X/2]]
continuous_events_4 = [(A ~ t_4) => [A ~ A + 2.0, X ~ X/2]]
discrete_events_1 = [
    10.0 => [X2_1 ~ X2_1 + 1.0]
    [5.0, 10.0] => [b_1 ~ 2 * b_1]
    (X > 5.0) => [X2_1 ~ X2_1 + 1.0, X ~ X - 1]
]
discrete_events_2 = [
    10.0 => [X2_2 ~ X2_2 + 1.0]
    [5.0, 10.0] => [b_2 ~ 2 * b_2]
    (X > 5.0) => [X2_2 ~ X2_2 + 1.0, X ~ X - 1]
]
discrete_events_3 = [
    10.0 => [X2_3 ~ X2_3 + 1.0]
    [5.0, 10.0] => [b_3 ~ 2 * b_3]
    (X > 5.0) => [X2_3 ~ X2_3 + 1.0, X ~ X - 1]
]
discrete_events_4 = [
    10.0 => [X2_4 ~ X2_4 + 1.0]
    [5.0, 10.0] => [b_4 ~ 2 * b_4]
    (X > 5.0) => [X2_4 ~ X2_4 + 1.0, X ~ X - 1]
]

# Creates the systems.
@named rs_4 = ReactionSystem(eqs_4, t; observed = observed_4, continuous_events = continuous_events_4,
                            discrete_events = discrete_events_4, spatial_ivs = sivs, 
                            metadata = "System 4", systems = [])
@named rs_2 = ReactionSystem(eqs_2, t; observed = observed_2, continuous_events = continuous_events_2,
                            discrete_events = discrete_events_2, spatial_ivs = sivs, 
                            metadata = "System 2", systems = [])
@named rs_3 = ReactionSystem(eqs_3, t; observed = observed_3, continuous_events = continuous_events_3,
                            discrete_events = discrete_events_3, spatial_ivs = sivs, 
                            metadata = "System 3", systems = [rs_4])
@named rs_1 = ReactionSystem(eqs_1, t; observed = observed_1, continuous_events = continuous_events_1,
                            discrete_events = discrete_events_1, spatial_ivs = sivs, 
                            metadata = "System 1", systems = [rs_2, rs_3])
rs = complete(rs_1)

save_reaction_network("file.jl", rs)

Is printed to this file:

let

# Independent variable:
@variables t

# Spatial independent variables:
spatial_ivs = @variables x y z [description = "A spatial independent variable."]

# Parameters:
ps = @parameters p d k1_1 k2_1 a b_1 t_1=2.0

# Species:
sps = @species X(t) X2_1(t)

# Variables:
vars = @variables A(t)=p [description = "A variable."] B_1(t)

# Reactions:
rxs = [
	Reaction(p, nothing, [X], nothing, [1]; metadata = [:description => "A reaction"]),
	Reaction(d, [X], nothing, [1], nothing; metadata = [:noise_scaling => h]),
	Reaction(k1_1, [X], [X2_1], [2], [1]),
	Reaction(k2_1, [X2_1], [X], [1], [2])
]

# Equations:
eqs = [
	Differential(t)(A) ~ -A + a,
	A + 2(B_1^3) ~ X*b_1
]

# Observables:
@variables O_1(t)
observed = [O_1 ~ X + 2X2_1]

# Continuous events:
continuous_events = [[A ~ t_1] => [A ~ 2.0 + A, X ~ (1//2)*X]]

# Discrete events:
discrete_events = [
	10.0 => [X2_1 ~ 1.0 + X2_1],
	[5.0, 10.0] => [b_1 ~ 2b_1],
	(X > 5.0) => [X2_1 ~ 1.0 + X2_1, X ~ -1 + X]
]

# Subystems:
systems = Vector(undef, 2)

# Declares subsystem: rs_2
systems[1] = let
	# Independent variable:
	@variables t
	
	# Spatial independent variables:
	local spatial_ivs = @variables x y z [description = "A spatial independent variable."]
	
	# Parameters:
	local ps = @parameters p d k1_2 k2_2 a b_2 t_2::Float64
	
	# Species:
	local sps = @species X(t) X2_2(t)
	
	# Variables:
	local vars = @variables A(t)=p [description = "A variable."] B_2(t)
	
	# Reactions:
	local rxs = [
		Reaction(p, nothing, [X], nothing, [1]; metadata = [:description => "A reaction"]),
		Reaction(d, [X], nothing, [1], nothing; metadata = [:noise_scaling => h]),
		Reaction(k1_2, [X], [X2_2], [2], [1]),
		Reaction(k2_2, [X2_2], [X], [1], [2])
	]
	
	# Equations:
	local eqs = [
		Differential(t)(A) ~ -A + a,
		A + 2(B_2^3) ~ X*b_2
	]
	
	# Observables:
	@variables O_2(t)
	local observed = [O_2 ~ X + 2X2_2]
	
	# Continuous events:
	local continuous_events = [[A ~ t_2] => [A ~ 2.0 + A, X ~ (1//2)*X]]
	
	# Discrete events:
	local discrete_events = [
		10.0 => [X2_2 ~ 1.0 + X2_2],
		[5.0, 10.0] => [b_2 ~ 2b_2],
		(X > 5.0) => [X2_2 ~ 1.0 + X2_2, X ~ -1 + X]
	]
	
	# Declares ReactionSystem model:
	ReactionSystem([rxs; eqs], t, [sps; vars], ps; name = :rs_2, spatial_ivs, observed, continuous_events, discrete_events, metadata = "System 2")
end

# Declares subsystem: rs_3
systems[2] = let
	# Independent variable:
	@variables t
	
	# Spatial independent variables:
	local spatial_ivs = @variables x y z [description = "A spatial independent variable."]
	
	# Parameters:
	local ps = @parameters p d k1_3 k2_3 a b_3 t_3 [description = "A parameter."]
	
	# Species:
	local sps = @species X(t) X2_3(t)
	
	# Variables:
	local vars = @variables A(t)=p [description = "A variable."] B_3(t)
	
	# Reactions:
	local rxs = [
		Reaction(p, nothing, [X], nothing, [1]; metadata = [:description => "A reaction"]),
		Reaction(d, [X], nothing, [1], nothing; metadata = [:noise_scaling => h]),
		Reaction(k1_3, [X], [X2_3], [2], [1]),
		Reaction(k2_3, [X2_3], [X], [1], [2])
	]
	
	# Equations:
	local eqs = [
		Differential(t)(A) ~ -A + a,
		A + 2(B_3^3) ~ X*b_3
	]
	
	# Observables:
	@variables O_3(t)
	local observed = [O_3 ~ X + 2X2_3]
	
	# Continuous events:
	local continuous_events = [[A ~ t_3] => [A ~ 2.0 + A, X ~ (1//2)*X]]
	
	# Discrete events:
	local discrete_events = [
		10.0 => [X2_3 ~ 1.0 + X2_3],
		[5.0, 10.0] => [b_3 ~ 2b_3],
		(X > 5.0) => [X2_3 ~ 1.0 + X2_3, X ~ -1 + X]
	]
	
	# Subystems:
	local systems = Vector(undef, 1)
	
	# Declares subsystem: rs_4
	systems[1] = let
		# Independent variable:
		@variables t
		
		# Spatial independent variables:
		local spatial_ivs = @variables x y z [description = "A spatial independent variable."]
		
		# Species:
		local sps = @species X(t) X2_4(t)=p [description = "A species."]
		
		# Variables:
		local vars = @variables A(t)=p [description = "A variable."] B_4(t)
		
		# Some parameters depends on the declaration of other parameters, species, and/or variables.
		# These are specially handled here.
		@parameters p d k1_4 k2_4 a b_4
		@parameters t_4::Float32=p [description = "A parameter."]
		local ps = [p, d, k1_4, k2_4, a, b_4, t_4]
		
		# Reactions:
		local rxs = [
			Reaction(p, nothing, [X], nothing, [1]; metadata = [:description => "A reaction"]),
			Reaction(d, [X], nothing, [1], nothing; metadata = [:noise_scaling => h]),
			Reaction(k1_4, [X], [X2_4], [2], [1]),
			Reaction(k2_4, [X2_4], [X], [1], [2])
		]
		
		# Equations:
		local eqs = [
			Differential(t)(A) ~ -A + a,
			A + 2(B_4^3) ~ X*b_4
		]
		
		# Observables:
		@variables O_4(t)
		local observed = [O_4 ~ X + 2X2_4]
		
		# Continuous events:
		local continuous_events = [[A ~ t_4] => [A ~ 2.0 + A, X ~ (1//2)*X]]
		
		# Discrete events:
		local discrete_events = [
			10.0 => [X2_4 ~ 1.0 + X2_4],
			[5.0, 10.0] => [b_4 ~ 2b_4],
			(X > 5.0) => [X2_4 ~ 1.0 + X2_4, X ~ -1 + X]
		]
		
		# Declares ReactionSystem model:
		ReactionSystem([rxs; eqs], t, [sps; vars], ps; name = :rs_4, spatial_ivs, observed, continuous_events, discrete_events, metadata = "System 4")
	end
	
	# Declares ReactionSystem model:
	ReactionSystem([rxs; eqs], t, [sps; vars], ps; name = :rs_3, spatial_ivs, observed, continuous_events, discrete_events, systems, metadata = "System 3")
end

# Declares ReactionSystem model:
rs = ReactionSystem([rxs; eqs], t, [sps; vars], ps; name = :rs_1, spatial_ivs, observed, continuous_events, discrete_events, systems, metadata = "System 1")
complete(rs)

end

[isaacsas]

How does this interact with non-traceable user-registered functions? (Things like using user-written functions in rates via @register_symbolic?)

[TorkelE]

Stuff like

f(x) = log(x) - 5x * sin(7)
# Declares the model.
rs = @reaction_network begin
    (f(p),d), 0 <--> X
end

works, however, I am not sure if there is anyway to deal with e.g.

g(x) = exp(x) - 1x
@register_symbolic g(x)
# Declares the model.
rs = @reaction_network begin
    (p,g(d)), 0 <--> X
end

[isaacsas]

Can you given an informative error in that case then?

[TorkelE]

I will see if I can come up with something