Merge branch 'master' into revert_testfile

Author: TorkelE

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 = [:missing_docs])
+    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

@@ -18,21 +18,20 @@ pages = Any[
         "Model creation examples" => Any[
             "model_creation/examples/basic_CRN_library.md",
             "model_creation/examples/programmatic_generative_linear_pathway.md",
-            "model_creation/examples/hodgkin_huxley_equation.md",
-            #"model_creation/examples/smoluchowski_coagulation_equation.md"
+            "model_creation/examples/hodgkin_huxley_equation.md"            #"model_creation/examples/smoluchowski_coagulation_equation.md"
         ]
     ],
     "Model simulation" => Any[
         "model_simulation/simulation_introduction.md",
         "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",
         "steady_state_functionality/nonlinear_solve.md",
-        "steady_state_functionality/steady_state_stability_computation.md",        
+        "steady_state_functionality/steady_state_stability_computation.md",
         "steady_state_functionality/bifurcation_diagrams.md",
         "steady_state_functionality/dynamical_systems.md"
     ],

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

ext/CatalystStructuralIdentifiabilityExtension/structural_identifiability_extension.jl

@@ -26,12 +26,13 @@ Notes:
 - This function is part of the StructuralIdentifiability.jl extension. StructuralIdentifiability.jl must be imported to access it.
 - `measured_quantities` and `known_p` input may also be symbolic (e.g. measured_quantities = [rs.X])
 """
-function Catalyst.make_si_ode(rs::ReactionSystem; measured_quantities = [], known_p = [], 
-                              ignore_no_measured_warn = false, remove_conserved = true)
+function Catalyst.make_si_ode(rs::ReactionSystem; measured_quantities = [], known_p = [],
+        ignore_no_measured_warn = false, remove_conserved = true)
     # Creates a MTK ODESystem, and a list of measured quantities (there are equations).
     # Gives these to SI to create an SI ode model of its preferred form.
     osys, conseqs, _ = make_osys(rs; remove_conserved)
-    measured_quantities = make_measured_quantities(rs, measured_quantities, known_p, conseqs; ignore_no_measured_warn)
+    measured_quantities = make_measured_quantities(rs, measured_quantities, known_p,
+        conseqs; ignore_no_measured_warn)
     return SI.mtk_to_si(osys, measured_quantities)[1]
 end
 
@@ -62,16 +63,18 @@ Notes:
 - This function is part of the StructuralIdentifiability.jl extension. StructuralIdentifiability.jl must be imported to access it.
 - `measured_quantities` and `known_p` input may also be symbolic (e.g. measured_quantities = [rs.X])
 """
-function SI.assess_local_identifiability(rs::ReactionSystem, args...; measured_quantities = [], 
-                                         known_p = [], funcs_to_check = Vector(), remove_conserved = true, 
-                                         ignore_no_measured_warn=false, kwargs...)
+function SI.assess_local_identifiability(rs::ReactionSystem, args...;
+        measured_quantities = [], known_p = [], funcs_to_check = Vector(),
+        remove_conserved = true, ignore_no_measured_warn = false, kwargs...)
     # Creates a ODESystem, list of measured quantities, and functions to check, of SI's preferred form.
     osys, conseqs, vars = make_osys(rs; remove_conserved)
-    measured_quantities = make_measured_quantities(rs, measured_quantities, known_p, conseqs; ignore_no_measured_warn)
+    measured_quantities = make_measured_quantities(rs, measured_quantities, known_p,
+        conseqs; ignore_no_measured_warn)
     funcs_to_check = make_ftc(funcs_to_check, conseqs, vars)
 
     # Computes identifiability and converts it to a easy to read form.
-    out = SI.assess_local_identifiability(osys, args...; measured_quantities, funcs_to_check, kwargs...)
+    out = SI.assess_local_identifiability(sys, args...; measured_quantities,
+        funcs_to_check, kwargs...)
     return make_output(out, funcs_to_check, reverse.(conseqs))
 end
 
@@ -100,16 +103,18 @@ Notes:
 - This function is part of the StructuralIdentifiability.jl extension. StructuralIdentifiability.jl must be imported to access it.
 - `measured_quantities` and `known_p` input may also be symbolic (e.g. measured_quantities = [rs.X])
 """
-function SI.assess_identifiability(rs::ReactionSystem, args...; measured_quantities = [], known_p = [], 
-                                   funcs_to_check = Vector(), remove_conserved = true, 
-                                   ignore_no_measured_warn=false, kwargs...)
+function SI.assess_identifiability(rs::ReactionSystem, args...;
+        measured_quantities = [], known_p = [], funcs_to_check = Vector(),
+        remove_conserved = true, ignore_no_measured_warn = false, kwargs...)
     # Creates a ODESystem, list of measured quantities, and functions to check, of SI's preferred form.
     osys, conseqs, vars = make_osys(rs; remove_conserved)
-    measured_quantities = make_measured_quantities(rs, measured_quantities, known_p, conseqs; ignore_no_measured_warn)
+    measured_quantities = make_measured_quantities(rs, measured_quantities, known_p,
+        conseqs; ignore_no_measured_warn)
     funcs_to_check = make_ftc(funcs_to_check, conseqs, vars)
 
     # Computes identifiability and converts it to a easy to read form.
-    out = SI.assess_identifiability(osys, args...; measured_quantities, funcs_to_check, kwargs...)
+    out = SI.assess_identifiability(osys, args...; measured_quantities,
+        funcs_to_check, kwargs...)
     return make_output(out, funcs_to_check, reverse.(conseqs))
 end
 
@@ -138,12 +143,13 @@ Notes:
 - This function is part of the StructuralIdentifiability.jl extension. StructuralIdentifiability.jl must be imported to access it.
 - `measured_quantities` and `known_p` input may also be symbolic (e.g. measured_quantities = [rs.X])
 """
-function SI.find_identifiable_functions(rs::ReactionSystem, args...; measured_quantities = [], 
-                                        known_p = [], remove_conserved = true, ignore_no_measured_warn=false, 
-                                        kwargs...)
+function SI.find_identifiable_functions(rs::ReactionSystem, args...;
+        measured_quantities = [], known_p = [], remove_conserved = true,
+        ignore_no_measured_warn = false, kwargs...)
     # Creates a ODESystem, and list of measured quantities, of SI's preferred form.
     osys, conseqs = make_osys(rs; remove_conserved)
-    measured_quantities = make_measured_quantities(rs, measured_quantities, known_p, conseqs; ignore_no_measured_warn)
+    measured_quantities = make_measured_quantities(rs, measured_quantities, known_p,
+        conseqs; ignore_no_measured_warn)
 
     # Computes identifiable functions and converts it to a easy to read form.
     out = SI.find_identifiable_functions(osys, args...; measured_quantities, kwargs...)
@@ -154,10 +160,12 @@ end
 
 # From a reaction system, creates the corresponding MTK-style ODESystem for SI application 
 # Also compute the, later needed, conservation law equations and list of system symbols (unknowns and parameters).
-function make_osys(rs::ReactionSystem; remove_conserved=true)
+function make_osys(rs::ReactionSystem; remove_conserved = true)
     # Creates the ODESystem corresponding to the ReactionSystem (expanding functions and flattening it).
     # Creates a list of the systems all symbols (unknowns and parameters).
-    ModelingToolkit.iscomplete(rs) || error("Identifiability should only be computed for complete systems. A ReactionSystem can be marked as complete using the `complete` function.")
+    if !ModelingToolkit.iscomplete(rs)
+        error("Identifiability should only be computed for complete systems. A ReactionSystem can be marked as complete using the `complete` function.")
+    end
     rs = complete(Catalyst.expand_registered_functions(flatten(rs)))
     osys = complete(convert(ODESystem, rs; remove_conserved))
     vars = [unknowns(rs); parameters(rs)]
@@ -177,10 +185,11 @@ end
 # Creates a list of measured quantities of a form that SI can read.
 # Each measured quantity must have a form like:
 # `obs_var ~ X` # (Here, `obs_var` is a variable, and X is whatever we can measure).
-function make_measured_quantities(rs::ReactionSystem, measured_quantities::Vector{T}, known_p::Vector{S},
-                                  conseqs; ignore_no_measured_warn=false) where {T,S}
+function make_measured_quantities(
+        rs::ReactionSystem, measured_quantities::Vector{T}, known_p::Vector{S},
+        conseqs; ignore_no_measured_warn = false) where {T, S}
     # Warning if the user didn't give any measured quantities.
-    if ignore_no_measured_warn || isempty(measured_quantities) 
+    if ignore_no_measured_warn || isempty(measured_quantities)
         @warn "No measured quantity provided to the `measured_quantities` argument, any further identifiability analysis will likely fail. You can disable this warning by setting `ignore_no_measured_warn=true`."
     end
 
@@ -192,7 +201,8 @@ function make_measured_quantities(rs::ReactionSystem, measured_quantities::Vecto
     # Creates one internal observation variable for each measured quantity (`___internal_observables`).
     # Creates a vector of equations, setting each measured quantity equal to one observation variable.
     @variables t (___internal_observables(Catalyst.get_iv(rs)))[1:length(mqs)]
-    return Equation[(q isa Equation) ? q : (___internal_observables[i] ~ q) for (i,q) in enumerate(mqs)] 
+    return Equation[(q isa Equation) ? q : (___internal_observables[i] ~ q)
+                    for (i, q) in enumerate(mqs)]
 end
 
 # Creates the functions that we wish to check for identifiability.
@@ -211,9 +221,9 @@ end
 function make_output(out, funcs_to_check, conseqs)
     funcs_to_check = vector_subs(funcs_to_check, conseqs)
     out = Dict(zip(vector_subs(keys(out), conseqs), values(out)))
-    sortdict = Dict(ftc => i for (i,ftc) in enumerate(funcs_to_check))
-    return sort(out; by = x -> sortdict[x])    
+    sortdict = Dict(ftc => i for (i, ftc) in enumerate(funcs_to_check))
+    return sort(out; by = x -> sortdict[x])
 end
 
 # For a vector of expressions and a conservation law, substitutes the law into every equation.
-vector_subs(eqs, subs) = [substitute(eq, subs) for eq in eqs]
+vector_subs(eqs, subs) = [substitute(eq, subs) for eq in eqs]

src/Catalyst.jl

@@ -7,7 +7,7 @@ using DocStringExtensions
 using SparseArrays, DiffEqBase, Reexport, Setfield
 using LaTeXStrings, Latexify, Requires
 using LinearAlgebra, Combinatorics
-using JumpProcesses: JumpProcesses, JumpProblem, 
+using JumpProcesses: JumpProcesses, JumpProblem,
                      MassActionJump, ConstantRateJump, VariableRateJump,
                      SpatialMassActionJump
 
@@ -16,7 +16,6 @@ using ModelingToolkit
 const MT = ModelingToolkit
 using DynamicQuantities#, Unitful # Having Unitful here as well currently gives an error.
 
-
 @reexport using ModelingToolkit
 using Symbolics
 using LinearAlgebra
@@ -81,7 +80,7 @@ const CONSERVED_CONSTANT_SYMBOL = :Γ
 # Declares symbols which may neither be used as parameters nor unknowns.
 const forbidden_symbols_skip = Set([:ℯ, :pi, :π, :t, :∅])
 const forbidden_symbols_error = union(Set([:im, :nothing, CONSERVED_CONSTANT_SYMBOL]),
-                                      forbidden_symbols_skip)
+    forbidden_symbols_skip)
 const forbidden_variables_error = let
     fvars = copy(forbidden_symbols_error)
     delete!(fvars, :t)
@@ -183,7 +182,6 @@ include("spatial_reaction_systems/utility.jl")
 include("spatial_reaction_systems/spatial_ODE_systems.jl")
 include("spatial_reaction_systems/lattice_jump_systems.jl")
 
-
 ### ReactionSystem Serialisation ###
 # Has to be at the end (because it uses records of all metadata declared by Catalyst).
 include("reactionsystem_serialisation/serialisation_support.jl")