Catalyst.jl
Catalyst.jl copied to clipboard
Published
20 hours ago •
SciML
SciML
Serialise `ReactionSystem` to file.
trafficstars
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
How does this interact with non-traceable user-registered functions? (Things like using user-written functions in rates via @register_symbolic?)
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