Skip to content

Commit

Permalink
updates
Browse files Browse the repository at this point in the history
  • Loading branch information
isaacsas committed Mar 2, 2023
1 parent 9aee17c commit b54e32a
Show file tree
Hide file tree
Showing 7 changed files with 75 additions and 10 deletions.
21 changes: 21 additions & 0 deletions HISTORY.md
Original file line number Diff line number Diff line change
Expand Up @@ -2,6 +2,27 @@

## Catalyst unreleased (master branch)

## Catalyst 13.1
- Non-species states can be declared in the DSL using `@variables`, and custom
independent variables (instead of just `t`) using `@ivs`. For the latter, the
first independent variable is always interpreted as the time variable, and all
*discovered* species are created to be functions of all the `ivs`. For example
in
```julia
rn = @reaction_network begin
@ivs s x
@variables A(s) B(x) C(s,x)
@species D(s) E(x) F(s,x)
k*C, A*D + B*E --> F + H
end
```
`s` will be the time variable, `H = H(s,x)` will be made a function of `s` and
`x`, and `A(s)`, `B(x)`, and `C(s,x)` will be non-species state variables.
- `Catalyst.isequal_ignore_names` has been deprecated for `isequivalent(rn1,
rn2)` to test equality of two networks and ignore their name. To include names
in the equality check continue to use `rn1 == rn2` or use `isequivalent(rn1,
rn2; ignorenames = false)`.

## Catalyst 13.0
- **BREAKING:** Parameters should no longer be listed at the end of the DSL
macro, but are instead inferred from their position in the reaction statements
Expand Down
2 changes: 1 addition & 1 deletion docs/src/api/catalyst_api.md
Original file line number Diff line number Diff line change
Expand Up @@ -228,7 +228,7 @@ reset_networkproperties!
## Network comparison
```@docs
==(rn1::Reaction, rn2::Reaction)
isequal_ignore_names
isequivalent
==(rn1::ReactionSystem, rn2::ReactionSystem)
```

Expand Down
36 changes: 35 additions & 1 deletion docs/src/introduction_to_catalyst/introduction_to_catalyst.md
Original file line number Diff line number Diff line change
Expand Up @@ -236,14 +236,48 @@ sprob = SDEProblem(bdp, u₀, tspan, p)
# solve and plot, tstops is used to specify enough points
# that the plot looks well-resolved
sol = solve(sprob, LambaEM(), tstops=range(0., step=4e-3, length=1001))
sol = solve(sprob, LambaEM(), tstops = range(0., step = 4e-3, length = 1001))
plot(sol)
```

We again have complete freedom to select any of the
StochasticDiffEq.jl SDE solvers, see the
[documentation](https://docs.sciml.ai/stable/modules/DiffEqDocs/solvers/sde_solve/).

---
## Specifying a complete model via the DSL
In the previous examples we specified initial conditions and parameter values
via mappings that were constructed after building our [`ReactionSystem`](@ref).
Catalyst also supports specifying default values for these during
`ReactionSystem` construction. For example, for the last SDE example we
could have also built and simulated the complete model using the DSL like
```@example tut1
bdp2 = @reaction_network begin
@parameters c₁ = 1.0 c₂ = 2.0 c₃ = 50.0
@species X(t) = 5.0
c₁, X --> 2X
c₂, X --> 0
c₃, 0 --> X
end
tspan = (0., 4.)
sprob2 = SDEProblem(bdp2, [], tspan)
```
Let's now simulate both models, starting from the same random number generator
seed, and check we get the same solutions
```@example tut1
using Random
Random.seed!(1)
sol = solve(sprob, LambaEM(), tstops = range(0., step = 4e-3, length = 1001))
p1 = plot(sol)
Random.seed!(1)
sol2 = solve(sprob2, LambaEM(), tstops = range(0., step = 4e-3, length = 1001))
p2 = plot(sol2)
plot(p1, p2, layout = (2,1))
```

For details on what information can be specified via the DSL see the [The
Reaction DSL](@ref dsl_description) tutorial.

---
## Reaction rate laws used in simulations
In generating mathematical models from a [`ReactionSystem`](@ref), reaction
Expand Down
1 change: 1 addition & 0 deletions src/Catalyst.jl
Original file line number Diff line number Diff line change
Expand Up @@ -75,6 +75,7 @@ export numspecies, numreactions, numreactionparams, setdefaults!, symmap_to_varm
export make_empty_network, addspecies!, addparam!, addreaction!, reactionparamsmap
export dependants, dependents, substoichmat, prodstoichmat, netstoichmat
export conservationlaws, conservedquantities, conservedequations, conservationlaw_constants
export isequivalent

# depreciated functions to remove in future releases
export params, numparams
Expand Down
17 changes: 13 additions & 4 deletions src/networkapi.jl
Original file line number Diff line number Diff line change
Expand Up @@ -1315,7 +1315,7 @@ function hash(rx::Reaction, h::UInt)
end

"""
isequal_ignore_names(rn1::ReactionSystem, rn2::ReactionSystem)
isequivalent(rn1::ReactionSystem, rn2::ReactionSystem; ignorenames = true)
Tests whether the underlying species, parameters and reactions are the same in
the two [`ReactionSystem`](@ref)s. Ignores the names of the systems in testing
Expand All @@ -1326,7 +1326,11 @@ Notes:
considered different than `(A+1)^2`.
- Does not include `defaults` in determining equality.
"""
function isequal_ignore_names(rn1::ReactionSystem, rn2::ReactionSystem)
function isequivalent(rn1::ReactionSystem, rn2::ReactionSystem; ignorenames = true)
if !ignorenames
(nameof(rn1) == nameof(rn2)) || return false
end

(get_combinatoric_ratelaws(rn1) == get_combinatoric_ratelaws(rn2)) || return false
isequal(get_iv(rn1), get_iv(rn2)) || return false
issetequal(get_sivs(rn1), get_sivs(rn2)) || return false
Expand All @@ -1342,6 +1346,12 @@ function isequal_ignore_names(rn1::ReactionSystem, rn2::ReactionSystem)
true
end

function isequal_ignore_names(rn1, rn2)
Base.depwarn("Catalyst.isequal_ignore_names has been deprecated. Use isequivalent(rn1, rn2) instead.",
:isequal_ignore_names; force = true)
isequivalent(rn1, rn2)
end

"""
==(rn1::ReactionSystem, rn2::ReactionSystem)
Expand All @@ -1355,8 +1365,7 @@ Notes:
- Does not include `defaults` in determining equality.
"""
function (==)(rn1::ReactionSystem, rn2::ReactionSystem)
(nameof(rn1) == nameof(rn2)) || return false
isequal_ignore_names(rn1, rn2)
isequivalent(rn1, rn2; ignorenames = false)
end

######################## functions to extend a network ####################
Expand Down
2 changes: 1 addition & 1 deletion src/reaction_network.jl
Original file line number Diff line number Diff line change
Expand Up @@ -370,7 +370,7 @@ function make_reaction_system(ex::Expr; name = :(gensym(:ReactionSystem)))
(sum(length.([reaction_lines, option_lines])) != length(ex.args)) &&
error("@reaction_network input contain $(length(ex.args) - sum(length.([reaction_lines,option_lines]))) malformed lines.")
any(!in(opt_in, option_keys) for opt_in in keys(options)) &&
error("The following unsupprted options were used: $(filter(opt_in->!in(opt_in,option_keys), keys(options)))")
error("The following unsupported options were used: $(filter(opt_in->!in(opt_in,option_keys), keys(options)))")
forbidden_symbol_check(union(species, parameters))
forbidden_variable_check(variables)
unique_symbol_check(union(species, parameters, variables, ivs))
Expand Down
6 changes: 3 additions & 3 deletions test/reactionsystem_structure/reactionsystem.jl
Original file line number Diff line number Diff line change
Expand Up @@ -92,7 +92,7 @@ end
# Test equation only constructor.
let
@named rs2 = ReactionSystem(rxs, t)
@test Catalyst.isequal_ignore_names(rs, rs2)
@test Catalyst.isequivalent(rs, rs2)
end

# Test show.
Expand Down Expand Up @@ -137,7 +137,7 @@ end
### Check ODE, SDE, and Jump Functions ###

# Test by evaluating drift and diffusion terms.
# Don't ask me (Torkel) why the statement before/after is needed.
# Don't ask me (Torkel) why the statement before/after is needed.
t = 0.0
let
p = rand(length(k))
Expand Down Expand Up @@ -384,7 +384,7 @@ let
rxs = [Reaction(x * t * A * B + y, [A], nothing)]
@named rs1 = ReactionSystem(rxs, t, [A, B], [x, y])
@named rs2 = ReactionSystem(rxs, t)
@test Catalyst.isequal_ignore_names(rs1, rs2)
@test Catalyst.isequivalent(rs1, rs2)

@variables t
@species L(t), H(t)
Expand Down

0 comments on commit b54e32a

Please sign in to comment.