diff --git a/Project.toml b/Project.toml index b831249319..1ac25b5040 100644 --- a/Project.toml +++ b/Project.toml @@ -1,6 +1,6 @@ name = "Catalyst" uuid = "479239e8-5488-4da2-87a7-35f2df7eef83" -version = "15.0.8" +version = "16" [deps] Combinatorics = "861a8166-3701-5b0c-9a16-15d98fcdc6aa" @@ -62,7 +62,7 @@ LaTeXStrings = "1.3.0" Latexify = "0.16.6" MacroTools = "0.5.5" Makie = "0.22.1" -ModelingToolkit = "9.73" +ModelingToolkit = "10" NetworkLayout = "0.4.7" Parameters = "0.12" Reexport = "1.0" diff --git a/src/Catalyst.jl b/src/Catalyst.jl index a079816d87..cedebfe319 100644 --- a/src/Catalyst.jl +++ b/src/Catalyst.jl @@ -36,7 +36,7 @@ import ModelingToolkit: get_variables, namespace_expr, namespace_equation, get_v # internal but needed ModelingToolkit functions import ModelingToolkit: check_variables, - check_parameters, _iszero, _merge, check_units, + check_parameters, _iszero, check_units, get_unit, check_equations, iscomplete import Base: (==), hash, size, getindex, setindex, isless, Sort.defalg, length, show diff --git a/src/reactionsystem.jl b/src/reactionsystem.jl index ad7f107961..31d3043afb 100644 --- a/src/reactionsystem.jl +++ b/src/reactionsystem.jl @@ -276,7 +276,7 @@ Notes: units). Unit checking can be disabled by passing the keyword argument `checks=false`. """ struct ReactionSystem{V <: NetworkProperties} <: - MT.AbstractTimeDependentSystem + MT.AbstractSystem """The equations (reactions and algebraic/differential) defining the system.""" eqs::Vector{CatalystEqType} """The Reactions defining the system. """ @@ -356,7 +356,7 @@ struct ReactionSystem{V <: NetworkProperties} <: check_parameters(ps, iv) nonrx_eqs = Equation[eq for eq in eqs if eq isa Equation] !isempty(nonrx_eqs) && check_equations(nonrx_eqs, iv) - !isempty(cevs) && check_equations(equations(cevs), iv) + !isnothing(cevs) && !isempty(cevs) && check_equations(equations(cevs), iv) end if isempty(sivs) && (checks == true || (checks & MT.CheckUnits) > 0) @@ -480,14 +480,10 @@ function ReactionSystem(eqs, iv, unknowns, ps; networkproperties end - # Creates the continuous and discrete callbacks. - ccallbacks = MT.SymbolicContinuousCallbacks(continuous_events) - dcallbacks = MT.SymbolicDiscreteCallbacks(discrete_events) - ReactionSystem( eqs′, rxs, iv′, sivs′, unknowns′, spcs, ps′, var_to_name, observed, name, systems, defaults, connection_type, nps, combinatoric_ratelaws, - ccallbacks, dcallbacks, metadata; checks = checks) + continuous_events, discrete_events, metadata; checks = checks) end # Two-argument constructor (reactions/equations and time variable). diff --git a/src/reactionsystem_conversions.jl b/src/reactionsystem_conversions.jl index 84d6060391..e9739a06ca 100644 --- a/src/reactionsystem_conversions.jl +++ b/src/reactionsystem_conversions.jl @@ -386,7 +386,7 @@ function assemble_jumps(rs; combinatoric_ratelaws = true, physical_scales = noth error("Must have at least one reaction that will be represented as a jump when constructing a JumpSystem.") # note isvrjvec indicates which reactions are not constant rate jumps - # it may be that a given jump has isvrjvec[i] = true but has a physical + # it may be that a given jump has isvrjvec[i] = true but has a physical isvrjvec = classify_vrjs(rs, physcales) rxvars = [] @@ -505,9 +505,16 @@ COMPLETENESS_ERROR = "A ReactionSystem must be complete before it can be convert ### System Conversions ### +abstract type SystemType end + +struct ReactionRateSystem <: SystemType end +struct ChemicalLangevinSystem <: SystemType end +struct StochasticChemicalKineticSystem <: SystemType end +struct NonlinearSystem <: SystemType end + """ ```julia -Base.convert(::Type{<:ODESystem},rs::ReactionSystem) +Base.convert(::Type{<:System},rs::ReactionSystem) ``` Convert a [`ReactionSystem`](@ref) to an `ModelingToolkit.ODESystem`. @@ -524,7 +531,7 @@ Keyword args and default values: with their rational function representation when converting to another system type. Set to `false`` to disable. """ -function Base.convert(::Type{<:ODESystem}, rs::ReactionSystem; name = nameof(rs), +function Base.convert(::Type{<:ReactionRateSystem}, rs::ReactionSystem; name = nameof(rs), combinatoric_ratelaws = get_combinatoric_ratelaws(rs), include_zero_odes = true, remove_conserved = false, checks = false, default_u0 = Dict(), default_p = Dict(), @@ -541,7 +548,7 @@ function Base.convert(::Type{<:ODESystem}, rs::ReactionSystem; name = nameof(rs) include_zero_odes, expand_catalyst_funs) eqs, us, ps, obs, defs = addconstraints!(eqs, fullrs, ists, ispcs; remove_conserved) - ODESystem(eqs, get_iv(fullrs), us, ps; + System(eqs, get_iv(fullrs), us, ps; observed = obs, name, defaults = _merge(defaults, defs), @@ -622,7 +629,7 @@ function Base.convert(::Type{<:NonlinearSystem}, rs::ReactionSystem; name = name all_differentials_permitted || nonlinear_convert_differentials_check(rs) eqs = [remove_diffs(eq.lhs) ~ remove_diffs(eq.rhs) for eq in eqs] - NonlinearSystem(eqs, us, ps; + System(eqs, us, ps; name, observed = obs, initialization_eqs = initeqs, defaults = _merge(defaults, defs), @@ -679,7 +686,7 @@ Notes: with their rational function representation when converting to another system type. Set to `false`` to disable. """ -function Base.convert(::Type{<:SDESystem}, rs::ReactionSystem; +function Base.convert(::Type{<:ChemicalLangevinSystem}, rs::ReactionSystem; name = nameof(rs), combinatoric_ratelaws = get_combinatoric_ratelaws(rs), include_zero_odes = true, checks = false, remove_conserved = false, default_u0 = Dict(), default_p = Dict(), @@ -704,7 +711,7 @@ function Base.convert(::Type{<:SDESystem}, rs::ReactionSystem; @info "Boundary condition species detected. As constraint equations are not currently supported when converting to SDESystems, the resulting system will be undetermined. Consider using constant species instead." end - SDESystem(eqs, noiseeqs, get_iv(flatrs), us, ps; + System([eqs; noiseeqs], get_iv(flatrs), us, ps; observed = obs, name, defaults = _merge(defaults, defs), @@ -728,7 +735,7 @@ Merge physical scales for a set of reactions. function merge_physical_scales(rxs, physical_scales, default) scales = get_physical_scale.(rxs) - # override metadata attached scales + # override metadata attached scales if physical_scales !== nothing for (key, scale) in physical_scales scales[key] = scale @@ -769,13 +776,13 @@ Notes: `VariableRateJump` to save the solution before and/or after the jump occurs. Defaults to true for both. """ -function Base.convert(::Type{<:JumpSystem}, rs::ReactionSystem; name = nameof(rs), +function Base.convert(::Type{<:StochasticChemicalKineticSystem}, rs::ReactionSystem; name = nameof(rs), combinatoric_ratelaws = get_combinatoric_ratelaws(rs), remove_conserved = nothing, checks = false, default_u0 = Dict(), default_p = Dict(), defaults = _merge(Dict(default_u0), Dict(default_p)), expand_catalyst_funs = true, save_positions = (true, true), physical_scales = nothing, kwargs...) iscomplete(rs) || error(COMPLETENESS_ERROR) - spatial_convert_err(rs::ReactionSystem, JumpSystem) + spatial_convert_err(rs::ReactionSystem, StochasticChemicalKineticSystem) (remove_conserved !== nothing) && throw(ArgumentError("Catalyst does not support removing conserved species when converting to JumpSystems.")) @@ -794,7 +801,7 @@ function Base.convert(::Type{<:JumpSystem}, rs::ReactionSystem; name = nameof(rs physical_scales, save_positions) ists, ispcs = get_indep_sts(flatrs) - # handle coupled ODEs and BC species + # handle coupled ODEs and BC species if (PhysicalScale.ODE in unique_scales) || has_nonreactions(flatrs) odeeqs = assemble_drift(flatrs, ispcs; combinatoric_ratelaws, remove_conserved = false, physical_scales, include_zero_odes = true) @@ -811,7 +818,7 @@ function Base.convert(::Type{<:JumpSystem}, rs::ReactionSystem; name = nameof(rs defs = MT.defaults(flatrs) end - JumpSystem(eqs, get_iv(flatrs), us, ps; + System(eqs, get_iv(flatrs), us, ps; observed = obs, name, defaults = _merge(defaults, defs), @@ -850,8 +857,8 @@ end DiffEqBase.NonlinearProblem(rs::ReactionSystem, u0, p = DiffEqBase.NullParameters(), args...; name = nameof(rs), combinatoric_ratelaws = get_combinatoric_ratelaws(rs), - remove_conserved = false, checks = false, check_length = false, - structural_simplify = remove_conserved, all_differentials_permitted = false, + remove_conserved = false, checks = false, check_length = false, + structural_simplify = remove_conserved, all_differentials_permitted = false, kwargs...) ``` @@ -923,7 +930,7 @@ Inputs for a JumpProblem from a given `ReactionSystem`. # Fields $(FIELDS) """ -struct JumpInputs{S <: MT.JumpSystem, T <: SciMLBase.AbstractODEProblem} +struct JumpInputs{S <: MT.System, T <: SciMLBase.AbstractODEProblem} """The `JumpSystem` to define the problem over""" sys::S """The problem the JumpProblem should be defined over, for example DiscreteProblem""" @@ -936,8 +943,8 @@ JumpInputs(rs::ReactionSystem, u0, tspan, p = DiffEqBase.NullParameters; name = nameof(rs), combinatoric_ratelaws = get_combinatoric_ratelaws(rs), - checks = false, physical_scales = nothing, - expand_catalyst_funs = true, + checks = false, physical_scales = nothing, + expand_catalyst_funs = true, save_positions = (true, true), remake_warn = true, kwargs...) ``` @@ -988,7 +995,7 @@ function JumpInputs(rs::ReactionSystem, u0, tspan, p = DiffEqBase.NullParameters name = nameof(rs), combinatoric_ratelaws = get_combinatoric_ratelaws(rs), checks = false, physical_scales = nothing, expand_catalyst_funs = true, save_positions = (true, true), remake_warn = true, kwargs...) - jsys = complete(convert(JumpSystem, rs; name, combinatoric_ratelaws, checks, + jsys = complete(convert(StochasticChemicalKineticSystem, rs; name, combinatoric_ratelaws, checks, physical_scales, expand_catalyst_funs, save_positions)) if MT.has_variableratejumps(jsys) || MT.has_equations(jsys) || @@ -1025,7 +1032,7 @@ function DiffEqBase.DiscreteProblem(rs::ReactionSystem, u0, tspan::Tuple, Base.depwarn("DiscreteProblem(rn::ReactionSystem, ...) is deprecated and will be \ removed in Catalyst 16. Use JumpInputs(rn, ...) instead.", :DiscreteProblem) - jsys = convert(JumpSystem, rs; name, combinatoric_ratelaws, checks, + jsys = convert(StochasticChemicalKineticSystem, rs; name, combinatoric_ratelaws, checks, expand_catalyst_funs) jsys = complete(jsys) return DiscreteProblem(jsys, u0, tspan, p, args...; kwargs...) @@ -1040,7 +1047,7 @@ function JumpProcesses.JumpProblem(rs::ReactionSystem, prob::SciMLBase.AbstractD Base.depwarn("JumpProblem(rn::ReactionSystem, prob, ...) is \ deprecated and will be removed in Catalyst 16. Use \ JumpProblem(JumpInputs(rn, ...), ...) instead.", :JumpProblem) - jsys = convert(JumpSystem, rs; name, combinatoric_ratelaws, + jsys = convert(StochasticChemicalKineticSystem, rs; name, combinatoric_ratelaws, expand_catalyst_funs, checks) jsys = complete(jsys) return JumpProblem(jsys, prob, aggregator; kwargs...) diff --git a/src/reactionsystem_serialisation/serialisation_support.jl b/src/reactionsystem_serialisation/serialisation_support.jl index bb6e7305ed..aec9fbdebc 100644 --- a/src/reactionsystem_serialisation/serialisation_support.jl +++ b/src/reactionsystem_serialisation/serialisation_support.jl @@ -223,7 +223,6 @@ const RECOGNISED_METADATA = Dict([Catalyst.ParameterConstantSpecies => "isconsta ModelingToolkit.VariableBounds => "bounds" ModelingToolkit.VariableUnit => "unit" ModelingToolkit.VariableConnectType => "connect" - ModelingToolkit.VariableNoiseType => "noise" ModelingToolkit.VariableInput => "input" ModelingToolkit.VariableOutput => "output" ModelingToolkit.VariableIrreducible => "irreducible" diff --git a/src/spatial_reaction_systems/lattice_reaction_systems.jl b/src/spatial_reaction_systems/lattice_reaction_systems.jl index e95c902176..33902af516 100644 --- a/src/spatial_reaction_systems/lattice_reaction_systems.jl +++ b/src/spatial_reaction_systems/lattice_reaction_systems.jl @@ -57,7 +57,7 @@ continuous space systems with them is possible, but requires the user to determi (the lattice). Better support for continuous space models is a work in progress. - Catalyst contains extensive documentation on spatial modelling, which can be found [here](https://docs.sciml.ai/Catalyst/stable/spatial_modelling/lattice_reaction_systems/). """ -struct LatticeReactionSystem{Q, R, S, T} <: MT.AbstractTimeDependentSystem +struct LatticeReactionSystem{Q, R, S, T} <: MT.AbstractSystem # Input values. """The (non-spatial) reaction system within each vertex.""" reactionsystem::ReactionSystem{Q} diff --git a/src/spatial_reaction_systems/utility.jl b/src/spatial_reaction_systems/utility.jl index 06e7e414d3..ceaf82c698 100644 --- a/src/spatial_reaction_systems/utility.jl +++ b/src/spatial_reaction_systems/utility.jl @@ -40,7 +40,7 @@ function lattice_process_p(ps_in, ps_vertex_syms::Vector, # For uniform parameters these have size 1/(1,1). Else, they have size num_verts/(num_verts,num_verts). ps = lattice_process_input(ps_in, [ps_vertex_syms; ps_edge_syms]) - # Split the parameter vector into one for vertex parameters and one for edge parameters. + # Split the parameter vector into one for vertex parameters and one for edge parameters. # Next, convert their values to the correct form (vectors for vert_ps and sparse matrices for edge_ps). vert_ps, edge_ps = split_parameters(ps, ps_vertex_syms, ps_edge_syms) vert_ps = vertex_value_map(vert_ps, lrs) @@ -164,13 +164,13 @@ function edge_value_form(values, lrs::LatticeReactionSystem, sym) throw(ArgumentError("Values was not provided for some edges for edge parameter $sym.")) end - # Unlike initial conditions/vertex parameters, (unless uniform) edge parameters' values are + # Unlike initial conditions/vertex parameters, (unless uniform) edge parameters' values are # always provided in the same (sparse matrix) form. return values end # Creates a map, taking each species (with transportation) to its transportation rate. -# The species is represented by its index (in species(lrs). +# The species is represented by its index (in species(lrs). # If the rate is uniform across all edges, the transportation rate will be a size (1,1) sparse matrix. # Else, the rate will be a size (num_verts,num_verts) sparse matrix. function make_sidxs_to_transrate_map(vert_ps::Vector{Pair{R, Vector{T}}}, @@ -190,7 +190,7 @@ end # Computes the transport rates for all species with transportation rates. Output is a map # taking each species' symbolics form to its transportation rates across all edges. function compute_all_transport_rates(p_val_dict, lrs::LatticeReactionSystem) - # For all species with transportation, compute their transportation rate (across all edges). + # For all species with transportation, compute their transportation rate (across all edges). # This is a vector, pairing each species to these rates. unsorted_rates = [s => compute_transport_rates(s, p_val_dict, lrs) for s in spatial_species(lrs)] @@ -199,7 +199,7 @@ function compute_all_transport_rates(p_val_dict, lrs::LatticeReactionSystem) return sort(unsorted_rates; by = rate -> findfirst(isequal(rate[1]), species(lrs))) end -# For the expression describing the rate of transport (likely only a single parameter, e.g. `D`), +# For the expression describing the rate of transport (likely only a single parameter, e.g. `D`), # and the values of all our parameters, compute the transport rate(s). # If all parameters that the rate depends on are uniform across all edges, this becomes a length-1 vector. # Else it becomes a vector where each value corresponds to the rate at one specific edge. @@ -304,10 +304,21 @@ end ### System Property Checks ### -# For a Symbolic expression, and a parameter set, check if any relevant parameters have a +# For a Symbolic expression, and a parameter set, check if any relevant parameters have a # spatial component. Filters out any parameters that are edge parameters. function has_spatial_vertex_component(exp, ps) relevant_syms = Symbolics.get_variables(exp) value_dict = Dict(filter(p -> p[2] isa Vector, ps)) return any(length(value_dict[sym]) > 1 for sym in relevant_syms) end + +### Utilities previously in ModelingToolkit.jl ### + +function todict(d) + eltype(d) <: Pair || throw(ArgumentError("The variable-value mapping must be a Dict.")) + d isa Dict ? d : Dict(d) +end + +_merge(d1, d2) = merge(todict(d1), todict(d2)) + +MT.refreshed_metadata(::Nothing) = MT.MetadataT() # FIXME: Type piracy