Skip to content

Commit

Permalink
Merge pull request #1006 from isaacsas/long_jump_plotting_fix
Browse files Browse the repository at this point in the history
Long jump plotting fix
  • Loading branch information
isaacsas committed Jul 25, 2024
2 parents cc1b958 + 0f8624f commit f99b6a6
Show file tree
Hide file tree
Showing 6 changed files with 113 additions and 40 deletions.
25 changes: 13 additions & 12 deletions docs/src/api.md
Original file line number Diff line number Diff line change
Expand Up @@ -68,11 +68,12 @@ jumpsys = convert(JumpSystem, rs)
jumpsys = complete(jumpsys)
u₀map = [S => 999, I => 1, R => 0]
dprob = DiscreteProblem(jumpsys, u₀map, tspan, parammap)
jprob = JumpProblem(jumpsys, dprob, Direct(); save_positions = (false,false))
sol = solve(jprob, SSAStepper(), saveat = 2.0)
jprob = JumpProblem(jumpsys, dprob, Direct())
sol = solve(jprob, SSAStepper())
sol = solve(jprob, SSAStepper(), seed = 1234) # hide
p3 = plot(sol, title = "jump")
plot(p1, p2, p3; layout = (3,1))
Catalyst.PNG(plot(p1, p2, p3; layout = (3,1), fmt = :png, dpi = 200)) # hide
```

```@docs
Expand All @@ -85,19 +86,19 @@ ReactionSystem
```

## [Options for the `@reaction_network` DSL](@id api_dsl_options)
We have [previously described](@ref dsl_advanced_options) how options permits the user to supply non-reaction information to [`ReactionSystem`](@ref) created through the DSL. Here follows a list
We have [previously described](@ref dsl_advanced_options) how options permits the user to supply non-reaction information to [`ReactionSystem`](@ref) created through the DSL. Here follows a list
of all options currently available.
- [`parameters`]:(@ref dsl_advanced_options_declaring_species_and_parameters) Allows the designation of a set of symbols as system parameter.
- [`parameters`]:(@ref dsl_advanced_options_declaring_species_and_parameters) Allows the designation of a set of symbols as system parameter.
- [`species`](@ref dsl_advanced_options_declaring_species_and_parameters): Allows the designation of a set of symbols as system species.
- [`variables`](@ref dsl_advanced_options_declaring_species_and_parameters): Allows the designation of a set of symbols as system non-species variables.
- [`variables`](@ref dsl_advanced_options_declaring_species_and_parameters): Allows the designation of a set of symbols as system non-species variables.
- [`ivs`](@ref dsl_advanced_options_ivs): Allows the designation of a set of symbols as system independent variables.
- [`compounds`](@ref chemistry_functionality_compounds): Allows the designation of compound species.
- [`observables`](@ref dsl_advanced_options_observables): Allows the designation of compound observables.
- [`observables`](@ref dsl_advanced_options_observables): Allows the designation of compound observables.
- [`default_noise_scaling`](@ref simulation_intro_SDEs_noise_saling): Enables the setting of a default noise scaling expression.
- [`differentials`](@ref constraint_equations_coupling_constraints): Allows the designation of differentials.
- [`equations`](@ref constraint_equations_coupling_constraints): Allows the creation of algebraic and/or differential equations.
- [`continuous_events`](@ref constraint_equations_events): Allows the creation of continuous events.
- [`discrete_events`](@ref constraint_equations_events): Allows the creation of discrete events.
- [`differentials`](@ref constraint_equations_coupling_constraints): Allows the designation of differentials.
- [`equations`](@ref constraint_equations_coupling_constraints): Allows the creation of algebraic and/or differential equations.
- [`continuous_events`](@ref constraint_equations_events): Allows the creation of continuous events.
- [`discrete_events`](@ref constraint_equations_events): Allows the creation of discrete events.
- [`combinatoric_ratelaws`](@ref faq_combinatoric_ratelaws): Takes a single option (`true` or `false`), which sets whether to use combinatorial rate laws.

## [ModelingToolkit and Catalyst accessor functions](@id api_accessor_functions)
Expand Down Expand Up @@ -186,7 +187,7 @@ isautonomous
```

## Coupled reaction/equation system properties
The following system property accessor functions are primarily relevant to reaction system [coupled
The following system property accessor functions are primarily relevant to reaction system [coupled
to differential and/or algebraic equations](@ref constraint_equations).
```@docs
ModelingToolkit.has_alg_equations
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -211,7 +211,7 @@ jprob = JumpProblem(rn, dprob, Direct())
# now, let's solve and plot the jump process:
sol = solve(jprob, SSAStepper())
plot(sol)
plot(sol, density = 10000, fmt = :png) # hide
Catalyst.PNG(plot(sol; fmt = :png, dpi = 200)) # hide
```

We see that oscillations remain, but become much noisier. Note, in constructing
Expand Down
48 changes: 23 additions & 25 deletions docs/src/model_creation/examples/basic_CRN_library.md
Original file line number Diff line number Diff line change
Expand Up @@ -116,11 +116,9 @@ oplt = plot(osol; title = "Reaction rate equation (ODE)")
splt = plot(ssol; title = "Chemical Langevin equation (SDE)")
jplt = plot(jsol; title = "Stochastic chemical kinetics (Jump)")
plot(oplt, splt, jplt; lw = 2, size=(800,800), layout = (3,1))
oplt = plot(osol; title = "Reaction rate equation (ODE)", plotdensity = 1000, fmt = :png) # hide
splt = plot(ssol; title = "Chemical Langevin equation (SDE)", plotdensity = 1000, fmt = :png) # hide
jplt = plot(jsol; title = "Stochastic chemical kinetics (Jump)", plotdensity = 1000, fmt = :png) # hide
plot(oplt, splt, jplt; lw = 2, size=(800,800), layout = (3,1), plotdensity = 1000, fmt = :png) # hide
plot!(bottom_margin = 3Plots.Measures.mm) # hide
fullplt = plot(oplt, splt, jplt; lw = 2, size=(800,800), layout = (3,1), fmt = :png,
dpi = 200, bottom_margin = 3Plots.Measures.mm) # hide
Catalyst.PNG(fullplt) # hide
```

## [SIR infection model](@id basic_CRN_library_sir)
Expand Down Expand Up @@ -148,23 +146,21 @@ Next, we perform 3 different Jump simulations. Note that for the stochastic mode
```@example crn_library_sir
using JumpProcesses
dprob = DiscreteProblem(sir_model, u0, tspan, ps)
jprob = JumpProblem(sir_model, dprob, Direct(); save_positions = (false, false))
jprob = JumpProblem(sir_model, dprob, Direct())
jsol1 = solve(jprob, SSAStepper(); saveat = 1.0)
jsol2 = solve(jprob, SSAStepper(); saveat = 1.0)
jsol3 = solve(jprob, SSAStepper(); saveat = 1.0)
jsol1 = solve(jprob, SSAStepper(); saveat = 1.0, seed = 1) # hide
jsol2 = solve(jprob, SSAStepper(); saveat = 1.0, seed = 2) # hide
jsol3 = solve(jprob, SSAStepper(); saveat = 1.0, seed = 3) # hide
jsol1 = solve(jprob, SSAStepper())
jsol2 = solve(jprob, SSAStepper())
jsol3 = solve(jprob, SSAStepper())
jsol1 = solve(jprob, SSAStepper(), seed = 1) # hide
jsol2 = solve(jprob, SSAStepper(), seed = 2) # hide
jsol3 = solve(jprob, SSAStepper(), seed = 3) # hide
jplt1 = plot(jsol1; title = "Outbreak")
jplt2 = plot(jsol2; title = "Outbreak")
jplt3 = plot(jsol3; title = "No outbreak")
plot(jplt1, jplt2, jplt3; lw = 3, size=(800,700), layout = (3,1))
jplt1 = plot(jsol1; title = "Outbreak", plotdensity = 1000, fmt = :png) # hide
jplt2 = plot(jsol2; title = "Outbreak", plotdensity = 1000, fmt = :png) # hide
jplt3 = plot(jsol3; title = "No outbreak", plotdensity = 1000, fmt = :png) # hide
plot(jplt1, jplt2, jplt3; lw = 3, size=(800,700), layout = (3,1), plotdensity = 1000, fmt = :png) # hide
fullplt = plot(jplt1, jplt2, jplt3; lw = 3, size=(800,700), layout = (3,1), fmt = :png, dpi = 200) # hide
Catalyst.PNG(fullplt) # hide
```

## [Chemical cross-coupling](@id basic_CRN_library_cc)
Expand All @@ -186,7 +182,7 @@ In two separate plots.
using OrdinaryDiffEq, Plots
u0 = [:S₁ => 1.0, :C => 0.05, :S₂ => 1.2, :S₁C => 0.0, :CP => 0.0, :P => 0.0]
tspan = (0., 15.)
ps = [:k₁ => 5.0, :k₂ => 5.0, :k₃ => 100.0]
ps = [:k₁ => 5.0, :k₂ => 5.0, :k₃ => 100.0]
# solve ODEs
oprob = ODEProblem(cc_system, u0, tspan, ps)
Expand All @@ -208,7 +204,7 @@ wilhelm_model = @reaction_network begin
k4, X --> 0
end
```
We can simulate the model for two different initial conditions, demonstrating the existence of two different stable steady states.
We can simulate the model for two different initial conditions, demonstrating the existence of two different stable steady states.
```@example crn_library_wilhelm
using OrdinaryDiffEq, Plots
u0_1 = [:X => 1.5, :Y => 0.5]
Expand All @@ -234,7 +230,7 @@ sa_loop = @reaction_network begin
d, X --> ∅
end
```
A simple example of such a loop is a transcription factor which activates its own gene. Here, $v₀$ represents a basic transcription rate (leakage) in the absence of the transcription factor.
A simple example of such a loop is a transcription factor which activates its own gene. Here, $v₀$ represents a basic transcription rate (leakage) in the absence of the transcription factor.

We simulate the self-activation loop from a single initial condition using both deterministic (ODE) and stochastic (jump) simulations. We note that while the deterministic simulation reaches a single steady state, the stochastic one switches between two different states.
```@example crn_library_self_activation
Expand All @@ -247,13 +243,15 @@ oprob = ODEProblem(sa_loop, u0, tspan, ps)
osol = solve(oprob)
dprob = DiscreteProblem(sa_loop, u0, tspan, ps)
jprob = JumpProblem(sa_loop, dprob, Direct(); save_positions = (false,false))
jsol = solve(jprob, SSAStepper(); saveat = 10.0)
jsol = solve(jprob, SSAStepper(); saveat = 10.0, seed = 12) # hide
jprob = JumpProblem(sa_loop, dprob, Direct())
jsol = solve(jprob, SSAStepper())
jsol = solve(jprob, SSAStepper(), seed = 12) # hide
plot(osol; lw = 3, label = "Reaction rate equation (ODE)")
plot!(jsol; lw = 3, label = "Stochastic chemical kinetics (Jump)", yguide = "X", size = (800,350))
plot!(bottom_margin = 3Plots.Measures.mm, left_margin=3Plots.Measures.mm) # hide
fplt = plot(osol; lw = 3, label = "Reaction rate equation (ODE)")
plot!(fplt, jsol; lw = 3, label = "Stochastic chemical kinetics (Jump)", yguide = "X", size = (800,350))
plot!(fplt, bottom_margin = 3Plots.Measures.mm, left_margin=3Plots.Measures.mm) # hide
plot!(fplt; fmt = :png, dpi = 200) # hide
Catalyst.PNG(fplt) # hide
```

## [The Brusselator](@id basic_CRN_library_brusselator)
Expand Down
2 changes: 1 addition & 1 deletion docs/src/model_creation/parametric_stoichiometry.md
Original file line number Diff line number Diff line change
Expand Up @@ -165,7 +165,7 @@ show(stdout, MIME"text/plain"(), equations(jsys)[4].rate) # hide
equations(jsys)[4].affect!
show(stdout, MIME"text/plain"(), equations(jsys)[4].affect!) # hide
```
Finally, we can now simulate our jumpsystem
Finally, we can now simulate our `JumpSystem`
```@example s1
pmean = 200
bval = 70
Expand Down
3 changes: 2 additions & 1 deletion src/Catalyst.jl
Original file line number Diff line number Diff line change
Expand Up @@ -141,7 +141,8 @@ export mm, mmr, hill, hillr, hillar
# for Latex printing of ReactionSystems
include("latexify_recipes.jl")

# for making and saving graphs
# for making and saving graphs/plots
include("plotting.jl")
include("graphs.jl")
export Graph, savegraph, complexgraph

Expand Down
73 changes: 73 additions & 0 deletions src/plotting.jl
Original file line number Diff line number Diff line change
@@ -0,0 +1,73 @@
#######################################################################################
"""
The following code is taken from the MIT licensed https://github.com/tkf/DisplayAs.jl
Copyright (c) 2019 Takafumi Arakaki <aka.tkf@gmail.com>
Permission is hereby granted, free of charge, to any person obtaining a copy of this
software and associated documentation files (the "Software"), to deal in the Software
without restriction, including without limitation the rights to use, copy, modify, merge,
publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons
to whom the Software is furnished to do so, subject to the following conditions:
The above copyright notice and this permission notice shall be included in all copies or
substantial portions of the Software.
THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED,
INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR
PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE
FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR
OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
DEALINGS IN THE SOFTWARE.
------------------------------------------------------------------------------------
In Catalyst this code is considered **internal**, and solely used for forcing a given image
type in documentation due to limitations in Documenter.jl and Plots.jl. It should never be
used for any other purpose or exported.
"""

# _showables = [
# (:Text, "text/plain")
# (:MD, "text/markdown")
# (:HTML, "text/html")
# (:JSON, "application/json")
# (:SVG, "image/svg+xml")
# (:PNG, "image/png")
# (:GIF, "image/gif")
# (:PDF, "application/pdf")
# (:EPS, "application/eps")
# (:JPEG, "image/jpeg")
# (:PS, "application/postscript")
# (:LaTeX, "text/latex")
# (:CSV, "text/csv")
# (:TSV, "text/tab-separated-values")
# ]

# currently we only use PNG conversion, but others can be added from above as needed
_showables = [(:PNG, "image/png")]

struct Showable{mime <: MIME}
content::Any
options::NamedTuple
end

function Showable{T}(content; options...) where {T <: MIME}
@nospecialize
return Showable{T}(content, (; options...))
end

for (_, mime) in _showables
MIMEType = typeof(MIME(mime))
@eval Base.show(io::IO, ::$MIMEType, s::Showable{>:$MIMEType}; options...) = show(
io, $MIMEType(), s.content; s.options..., options...)
end

macro mime_str(s)
:(Showable{MIME{Symbol($s)}})
end

for (name, mime) in _showables
@eval const $name = @mime_str $mime
end

#######################################################################################

0 comments on commit f99b6a6

Please sign in to comment.