merge master

Author: isaacsas

docs/make.jl

@@ -31,17 +31,17 @@ include("pages.jl")
 #          pages = pages)
 
 makedocs(sitename = "Catalyst.jl",
-         authors = "Samuel Isaacson",
-         format = Documenter.HTML(; analytics = "UA-90474609-3",
-                                  prettyurls = (get(ENV, "CI", nothing) == "true"),
-                                  assets = ["assets/favicon.ico"],
-                                  canonical = "https://docs.sciml.ai/Catalyst/stable/"),
-         modules = [Catalyst, ModelingToolkit],
-         doctest = false,
-         clean = true,
-         pages = pages,
-         pagesonly = true,
-         warnonly = true)
+    authors = "Samuel Isaacson",
+    format = Documenter.HTML(; analytics = "UA-90474609-3",
+        prettyurls = (get(ENV, "CI", nothing) == "true"),
+        assets = ["assets/favicon.ico"],
+        canonical = "https://docs.sciml.ai/Catalyst/stable/"),
+    modules = [Catalyst, ModelingToolkit],
+    doctest = false,
+    clean = true,
+    pages = pages,
+    pagesonly = true,
+    warnonly = [:missing_docs])
 
 deploydocs(repo = "github.com/SciML/Catalyst.jl.git";
-           push_preview = true)
+    push_preview = true)

docs/pages.jl

@@ -27,7 +27,7 @@ pages = Any[
         "model_simulation/simulation_plotting.md",
         "model_simulation/simulation_structure_interfacing.md",
         "model_simulation/ensemble_simulations.md",
-        "model_simulation/ode_simulation_performance.md",
+        "model_simulation/ode_simulation_performance.md"
     ],
     "Steady state analysis" => Any[
         "steady_state_functionality/homotopy_continuation.md",

docs/src/introduction_to_catalyst/introduction_to_catalyst.md

@@ -12,9 +12,9 @@ We first import the basic packages we'll need:
 
 ```@example tut1
 # If not already installed, first hit "]" within a Julia REPL. Then type:
-# add Catalyst DifferentialEquations Plots Latexify
+# add Catalyst OrdinaryDiffEq Plots Latexify
 
-using Catalyst, DifferentialEquations, Plots, Latexify
+using Catalyst, OrdinaryDiffEq, Plots, Latexify
 ```
 
 We now construct the reaction network. The basic types of arrows and predefined
@@ -160,7 +160,7 @@ underlying problem.
 
 At this point we are all set to solve the ODEs. We can now use any ODE solver
 from within the
-[DifferentialEquations.jl](https://docs.sciml.ai/DiffEqDocs/stable/solvers/ode_solve/)
+[OrdinaryDiffEq.jl](https://docs.sciml.ai/DiffEqDocs/stable/solvers/ode_solve/)
 package. We'll use the recommended default explicit solver, `Tsit5()`, and then
 plot the solutions:
 
@@ -169,7 +169,7 @@ sol = solve(oprob, Tsit5(), saveat=10.)
 plot(sol)
 ```
 We see the well-known oscillatory behavior of the repressilator! For more on the
-choices of ODE solvers, see the [DifferentialEquations.jl
+choices of ODE solvers, see the [OrdinaryDiffEq.jl
 documentation](https://docs.sciml.ai/DiffEqDocs/stable/solvers/ode_solve/).
 
 ---
@@ -182,6 +182,9 @@ Gillespie's `Direct` method, and then solve it to generate one realization of
 the jump process:
 
 ```@example tut1
+# imports the JumpProcesses packages 
+using JumpProcesses
+
 # redefine the initial condition to be integer valued
 u₀map = [:m₁ => 0, :m₂ => 0, :m₃ => 0, :P₁ => 20, :P₂ => 0, :P₃ => 0]
 
@@ -237,6 +240,9 @@ model by creating an `SDEProblem` and solving it similarly to what we did for OD
 above:
 
 ```@example tut1
+# imports the StochasticDiffEq package for SDE simulations
+using StochasticDiffEq
+
 # SDEProblem for CLE
 sprob = SDEProblem(bdp, u₀, tspan, p)
 

docs/src/model_creation/constraint_equations.md

@@ -28,7 +28,7 @@ creating these two systems.
 Here, to create differentials with respect to time (for our differential equations), we must import the time differential operator from Catalyst. We do this through `D = default_time_deriv()`. Here, `D(V)` denotes the differential of the variable `V` with respect to time.
 
 ```@example ceq1
-using Catalyst, DifferentialEquations, Plots
+using Catalyst, OrdinaryDiffEq, Plots
 t = default_t()
 D = default_time_deriv()
 
@@ -58,6 +58,7 @@ We can now merge the two systems into one complete `ReactionSystem` model using
 [`ModelingToolkit.extend`](@ref):
 ```@example ceq1
 @named growing_cell = extend(osys, rn)
+growing_cell = complete(growing_cell)
 ```
 
 We see that the combined model now has both the reactions and ODEs as its
@@ -72,7 +73,7 @@ plot(sol)
 As an alternative to the previous approach, we could have constructed our
 `ReactionSystem` all at once by directly using the symbolic interface:
 ```@example ceq2
-using Catalyst, DifferentialEquations, Plots
+using Catalyst, OrdinaryDiffEq, Plots
 t = default_t()
 D = default_time_deriv()
 
@@ -83,6 +84,7 @@ rx1 = @reaction $V, 0 --> P
 rx2 = @reaction 1.0, P --> 0
 @named growing_cell = ReactionSystem([rx1, rx2, eq], t)
 setdefaults!(growing_cell, [:P => 0.0])
+growing_cell = complete(growing_cell)
 
 oprob = ODEProblem(growing_cell, [], (0.0, 1.0))
 sol = solve(oprob, Tsit5())
@@ -100,12 +102,11 @@ the associated [ModelingToolkit
 tutorial](https://docs.sciml.ai/ModelingToolkit/stable/basics/Events/) for more
 details on the types of events that can be represented symbolically. A
 lower-level approach for creating events via the DifferentialEquations.jl
-callback interface is illustrated in the [Advanced Simulation Options](@ref
-advanced_simulations) tutorial.
+callback interface is illustrated [here](https://docs.sciml.ai/DiffEqDocs/stable/features/callback_functions/) tutorial.
 
 Let's first create our equations and unknowns/species again
 ```@example ceq3
-using Catalyst, DifferentialEquations, Plots
+using Catalyst, OrdinaryDiffEq, Plots
 t = default_t()
 D = default_time_deriv()
 
@@ -124,6 +125,7 @@ continuous_events = [V ~ 2.0] => [V ~ V/2, P ~ P/2]
 We can now create and simulate our model
 ```@example ceq3
 @named rs = ReactionSystem([rx1, rx2, eq], t; continuous_events)
+rs = complete(rs)
 
 oprob = ODEProblem(rs, [], (0.0, 10.0))
 sol = solve(oprob, Tsit5())
@@ -148,6 +150,7 @@ p = [k_on => 100.0, switch_time => 2.0, k_off => 10.0]
 Simulating our model, 
 ```@example ceq3
 @named osys = ReactionSystem(rxs, t, [A, B], [k_on, k_off, switch_time]; discrete_events)
+osys = complete(osys)
 
 oprob = ODEProblem(osys, u0, tspan, p)
 sol = solve(oprob, Tsit5(); tstops = 2.0)

ext/CatalystBifurcationKitExtension.jl

@@ -7,4 +7,4 @@ import BifurcationKit as BK
 # Creates and exports hc_steady_states function.
 include("CatalystBifurcationKitExtension/bifurcation_kit_extension.jl")
 
-end
+end

ext/CatalystBifurcationKitExtension/bifurcation_kit_extension.jl

@@ -2,23 +2,30 @@
 
 # Creates a BifurcationProblem, using a ReactionSystem as an input.
 function BK.BifurcationProblem(rs::ReactionSystem, u0_bif, ps, bif_par, args...;
-                               plot_var=nothing, record_from_solution=BK.record_sol_default, jac=true, u0=[], kwargs...)
-    if !isautonomous(rs) 
+        plot_var = nothing, record_from_solution = BK.record_sol_default, jac = true, u0 = [], kwargs...)
+    if !isautonomous(rs)
         error("Attempting to create a `BifurcationProblem` for a non-autonomous system (e.g. where some rate depend on $(rs.iv)). This is not possible.")
     end
 
     # Converts symbols to symbolics.
     (bif_par isa Symbol) && (bif_par = ModelingToolkit.get_var_to_name(rs)[bif_par])
     (plot_var isa Symbol) && (plot_var = ModelingToolkit.get_var_to_name(rs)[plot_var])
-    ((u0_bif isa Vector{<:Pair{Symbol,<:Any}}) || (u0_bif isa Dict{Symbol, <:Any})) && (u0_bif = symmap_to_varmap(rs, u0_bif))
-    ((ps isa Vector{<:Pair{Symbol,<:Any}}) || (ps isa Dict{Symbol, <:Any})) && (ps = symmap_to_varmap(rs, ps))
-    ((u0 isa Vector{<:Pair{Symbol,<:Any}}) || (u0 isa Dict{Symbol, <:Any})) && (u0 = symmap_to_varmap(rs, u0))
+    if (u0_bif isa Vector{<:Pair{Symbol, <:Any}}) || (u0_bif isa Dict{Symbol, <:Any})
+        u0_bif = symmap_to_varmap(rs, u0_bif)
+    end
+    if (ps isa Vector{<:Pair{Symbol, <:Any}}) || (ps isa Dict{Symbol, <:Any})
+        ps = symmap_to_varmap(rs, ps)
+    end
+    if (u0 isa Vector{<:Pair{Symbol, <:Any}}) || (u0 isa Dict{Symbol, <:Any})
+        u0 = symmap_to_varmap(rs, u0)
+    end
 
     # Creates NonlinearSystem.
     Catalyst.conservationlaw_errorcheck(rs, vcat(ps, u0))
-    nsys = complete(convert(NonlinearSystem, rs; remove_conserved=true, defaults=Dict(u0)))
+    nsys = convert(NonlinearSystem, rs; remove_conserved = true, defaults = Dict(u0))
+    nsys = complete(nsys)
 
     # Makes BifurcationProblem (this call goes through the ModelingToolkit-based BifurcationKit extension).
-    return BK.BifurcationProblem(nsys, u0_bif, ps, bif_par, args...; plot_var=plot_var,
-                                 record_from_solution=record_from_solution, jac=jac, kwargs...)
-end
+    return BK.BifurcationProblem(nsys, u0_bif, ps, bif_par, args...; plot_var,
+        record_from_solution, jac, kwargs...)
+end

ext/CatalystHomotopyContinuationExtension.jl

@@ -11,4 +11,4 @@ import Symbolics: unwrap, wrap, Rewriters, symtype, issym, istree
 # Creates and exports hc_steady_states function.
 include("CatalystHomotopyContinuationExtension/homotopy_continuation_extension.jl")
 
-end
+end

ext/CatalystHomotopyContinuationExtension/homotopy_continuation_extension.jl

@@ -35,8 +35,9 @@ Notes:
 - Homotopy-based steady state finding only works when all rates are rational polynomials (e.g. constant, linear, mm, or hill functions).
 ```
   """
-function Catalyst.hc_steady_states(rs::ReactionSystem, ps; filter_negative=true, neg_thres=-1e-20, u0=[], kwargs...)
-    if !isautonomous(rs) 
+function Catalyst.hc_steady_states(rs::ReactionSystem, ps; filter_negative = true,
+        neg_thres = -1e-20, u0 = [], kwargs...)
+    if !isautonomous(rs)
         error("Attempting to compute steady state for a non-autonomous system (e.g. where some rate depend on $(rs.iv)). This is not possible.")
     end
     ss_poly = steady_state_polynomial(rs, ps, u0)
@@ -49,10 +50,11 @@ end
 function steady_state_polynomial(rs::ReactionSystem, ps, u0)
     rs = Catalyst.expand_registered_functions(rs)
     ns = complete(convert(NonlinearSystem, rs; remove_conserved = true))
-    pre_varmap = [symmap_to_varmap(rs,u0)..., symmap_to_varmap(rs,ps)...]
+    pre_varmap = [symmap_to_varmap(rs, u0)..., symmap_to_varmap(rs, ps)...]
     Catalyst.conservationlaw_errorcheck(rs, pre_varmap)
-    p_vals = ModelingToolkit.varmap_to_vars(pre_varmap, parameters(ns); defaults = ModelingToolkit.defaults(ns))
-    p_dict  = Dict(parameters(ns) .=> p_vals)
+    p_vals = ModelingToolkit.varmap_to_vars(pre_varmap, parameters(ns);
+        defaults = ModelingToolkit.defaults(ns))
+    p_dict = Dict(parameters(ns) .=> p_vals)
     eqs_pars_funcs = vcat(equations(ns), conservedequations(rs))
     eqs = map(eq -> substitute(eq.rhs - eq.lhs, p_dict), eqs_pars_funcs)
     eqs_intexp = make_int_exps.(eqs)
@@ -61,12 +63,14 @@ function steady_state_polynomial(rs::ReactionSystem, ps, u0)
 end
 
 # Parses and expression and return a version where any exponents that are Float64 (but an int, like 2.0) are turned into Int64s.
-make_int_exps(expr) = wrap(Rewriters.Postwalk(Rewriters.PassThrough(___make_int_exps))(unwrap(expr))).val
+function make_int_exps(expr)
+    wrap(Rewriters.Postwalk(Rewriters.PassThrough(___make_int_exps))(unwrap(expr))).val
+end
 function ___make_int_exps(expr)
     !istree(expr) && return expr
-    if (operation(expr) == ^) 
+    if (operation(expr) == ^)
         if isinteger(arguments(expr)[2])
-            return arguments(expr)[1] ^ Int64(arguments(expr)[2])
+            return arguments(expr)[1]^Int64(arguments(expr)[2])
         else
             error("An non integer ($(arguments(expr)[2])) was found as a variable exponent. Non-integer exponents are not supported for homotopy continuation based steady state finding.")
         end
@@ -95,7 +99,7 @@ function reorder_sols!(sols, ss_poly, rs::ReactionSystem)
 end
 
 # Filters away solutions with negative species concentrations (and for neg_thres < val < 0.0, sets val=0.0).
-function filter_negative_f(sols; neg_thres=-1e-20)
+function filter_negative_f(sols; neg_thres = -1e-20)
     for sol in sols, idx in 1:length(sol)
         (neg_thres < sol[idx] < 0) && (sol[idx] = 0)
     end
@@ -104,9 +108,13 @@ end
 
 # Sometimes (when polynomials are created from coupled CRN/DAEs), the steady state polynomial have the wrong type.
 # This converts it to the correct type, which homotopy continuation can handle.
-const WRONG_POLY_TYPE = Vector{DynamicPolynomials.Polynomial{DynamicPolynomials.Commutative{DynamicPolynomials.CreationOrder}, DynamicPolynomials.Graded{DynamicPolynomials.LexOrder}}}
-const CORRECT_POLY_TYPE = Vector{DynamicPolynomials.Polynomial{DynamicPolynomials.Commutative{DynamicPolynomials.CreationOrder}, DynamicPolynomials.Graded{DynamicPolynomials.LexOrder}, Float64}}
+const WRONG_POLY_TYPE = Vector{DynamicPolynomials.Polynomial{
+    DynamicPolynomials.Commutative{DynamicPolynomials.CreationOrder},
+    DynamicPolynomials.Graded{DynamicPolynomials.LexOrder}}}
+const CORRECT_POLY_TYPE = Vector{DynamicPolynomials.Polynomial{
+    DynamicPolynomials.Commutative{DynamicPolynomials.CreationOrder},
+    DynamicPolynomials.Graded{DynamicPolynomials.LexOrder}, Float64}}
 function poly_type_convert(ss_poly)
     (typeof(ss_poly) == WRONG_POLY_TYPE) && return convert(CORRECT_POLY_TYPE, ss_poly)
     return ss_poly
-end
+end

ext/CatalystStructuralIdentifiabilityExtension.jl

@@ -7,4 +7,4 @@ import StructuralIdentifiability as SI
 # Creates and exports hc_steady_states function.
 include("CatalystStructuralIdentifiabilityExtension/structural_identifiability_extension.jl")
 
-end
+end