[WIP] v16 changes to get tests passing

.github/workflows/Test.yml

@@ -2,8 +2,6 @@ name: "Tests"
 
 on:
   pull_request:
-    branches:
-      - master
     paths-ignore:
       - 'docs/**'
   push:

Project.toml

@@ -62,7 +62,7 @@ LaTeXStrings = "1.3.0"
 Latexify = "0.16.6"
 MacroTools = "0.5.5"
 Makie = "0.22.1"
-ModelingToolkit = "9.73"
+ModelingToolkit = "9.73, 10"
 NetworkLayout = "0.4.7"
 Parameters = "0.12"
 Reexport = "1.0"

docs/src/api/core_api.md

@@ -64,7 +64,7 @@ sol = solve(sprob, EM(), dt=.01, saveat = 2.0)
 p2 = plot(sol, title = "SDE")
 
 # solve as jump process
-jumpsys = convert(JumpSystem, rs)
+jumpsys = make_sck_jump(rs)
 jumpsys = complete(jumpsys)
 u₀map    = [S => 999, I => 1, R => 0]
 dprob = DiscreteProblem(jumpsys, u₀map, tspan, parammap)

docs/src/model_creation/parametric_stoichiometry.md

@@ -148,7 +148,7 @@ The parameter `b` does not need to be explicitly declared in the
 We next convert our network to a jump process representation
 ```@example s1
 using JumpProcesses
-jsys = convert(JumpSystem, burstyrn; combinatoric_ratelaws = false)
+jsys = make_sck_jump(burstyrn; combinatoric_ratelaws = false)
 jsys = complete(jsys)
 equations(jsys)
 show(stdout, MIME"text/plain"(), equations(jsys)) # hide

ext/CatalystBifurcationKitExtension/bifurcation_kit_extension.jl

@@ -37,7 +37,6 @@ function bkext_make_nsys(rs, u0)
     cons_default = [cons_eq.rhs for cons_eq in cons_eqs]
     cons_default = Catalyst.get_networkproperties(rs).conservedconst => cons_default
     defaults = Dict([u0; cons_default])
-    nsys = convert(
-        NonlinearSystem, rs; defaults, remove_conserved = true, conseqs_remake_warn = false)
+    nsys = make_rre_algeqs(rs; defaults, remove_conserved = true, conseqs_remake_warn = false)
     return complete(nsys)
 end

ext/CatalystHomotopyContinuationExtension.jl

@@ -2,6 +2,7 @@ module CatalystHomotopyContinuationExtension
 
 # Fetch packages.
 using Catalyst
+import DiffEqBase
 import DynamicPolynomials
 import ModelingToolkit as MT
 import HomotopyContinuation as HC

ext/CatalystHomotopyContinuationExtension/homotopy_continuation_extension.jl

@@ -51,7 +51,7 @@ function steady_state_polynomial(rs::ReactionSystem, ps, u0)
     # Creates the appropriate nonlinear system, and converts parameters to a form that can
     # be substituted in later.
     rs = Catalyst.expand_registered_functions(rs)
-    ns = complete(convert(NonlinearSystem, rs; remove_conserved = true, conseqs_remake_warn = false))
+    ns = complete(make_rre_algeqs(rs; remove_conserved = true, conseqs_remake_warn = false))
     pre_varmap = [symmap_to_varmap(rs, u0)..., symmap_to_varmap(rs, ps)...]
     Catalyst.conservationlaw_errorcheck(rs, pre_varmap)
     p_dict = make_p_val_dict(pre_varmap, rs, ns)
@@ -82,7 +82,7 @@ function make_p_val_dict(pre_varmap, rs, ns)
     foreach(conseq -> defaults[conseq.lhs] = conseq.rhs, conservationlaw_constants(rs))
 
     # Creates and return the full parameter value dictionary.p_vals = ModelingToolkit.varmap_to_vars(pre_varmap, all_ps; defaults = def_dict)
-    p_vals = ModelingToolkit.varmap_to_vars(pre_varmap, ps; defaults)
+    p_vals = varmap_to_vars_mtkv9(pre_varmap, ps; defaults)
     return Dict(ps .=> p_vals)
 end
 
@@ -184,3 +184,95 @@ function poly_type_convert(ss_poly)
     (typeof(ss_poly) == WRONG_POLY_TYPE) && return convert(CORRECT_POLY_TYPE, ss_poly)
     return ss_poly
 end
+
+
+
+### SAVED ARCHIVED MTK FUNCTION - REMOVE SOME TIME ###
+# pre-v10 version of function
+function varmap_to_vars_mtkv9(varmap, varlist; defaults = Dict(), check = true,
+        toterm = ModelingToolkit.default_toterm, promotetoconcrete = nothing,
+        tofloat = true, use_union = true)
+    varlist = collect(map(unwrap, varlist))
+
+    # Edge cases where one of the arguments is effectively empty.
+    is_incomplete_initialization = varmap isa DiffEqBase.NullParameters ||
+                                   varmap === nothing
+    if is_incomplete_initialization || isempty(varmap)
+        if isempty(defaults)
+            if !is_incomplete_initialization && check
+                isempty(varlist) || throw(ModelingToolkit.MissingVariablesError(varlist))
+            end
+            return nothing
+        else
+            varmap = Dict()
+        end
+    end
+
+    # We respect the input type if it's a static array
+    # otherwise canonicalize to a normal array
+    # container_type = T <: Union{Dict,Tuple} ? Array : T
+    if varmap isa ModelingToolkit.StaticArray
+        container_type = typeof(varmap)
+    else
+        container_type = Array
+    end
+
+    vals = if eltype(varmap) <: Pair # `varmap` is a dict or an array of pairs
+        varmap = ModelingToolkit.todict(varmap)
+        _varmap_to_vars_mtkv9(varmap, varlist; defaults = defaults, check = check,
+            toterm = toterm)
+    else # plain array-like initialization
+        varmap
+    end
+
+    promotetoconcrete === nothing && (promotetoconcrete = container_type <: AbstractArray)
+    if promotetoconcrete
+        vals = ModelingToolkit.promote_to_concrete(vals; tofloat = tofloat, use_union = use_union)
+    end
+
+    if isempty(vals)
+        return nothing
+    elseif container_type <: Tuple
+        (vals...,)
+    else
+        SymbolicUtils.Code.create_array(container_type, eltype(vals), Val{1}(),
+            Val(length(vals)), vals...)
+    end
+end
+
+function _varmap_to_vars_mtkv9(varmap::Dict, varlist; defaults = Dict(), check = false,
+        toterm = Symbolics.diff2term, initialization_phase = false)
+    varmap = canonicalize_varmap_mtkv9(varmap; toterm)
+    defaults = canonicalize_varmap_mtkv9(defaults; toterm)
+    varmap = merge(defaults, varmap)
+    values = Dict()
+
+    T = Union{}
+    for var in varlist
+        var = ModelingToolkit.unwrap(var)
+        val = ModelingToolkit.unwrap(ModelingToolkit.fixpoint_sub(var, varmap; operator = Symbolics.Operator))
+        if !isequal(val, var)
+            values[var] = val
+        end
+    end
+    missingvars = setdiff(varlist, collect(keys(values)))
+    check && (isempty(missingvars) || throw(ModelingToolkit.MissingVariablesError(missingvars)))
+    return [values[ModelingToolkit.unwrap(var)] for var in varlist]
+end
+
+function canonicalize_varmap_mtkv9(varmap; toterm = Symbolics.diff2term)
+    new_varmap = Dict()
+    for (k, v) in varmap
+        k = ModelingToolkit.unwrap(k)
+        v = ModelingToolkit.unwrap(v)
+        new_varmap[k] = v
+        new_varmap[toterm(k)] = v
+        if Symbolics.isarraysymbolic(k) && Symbolics.shape(k) !== Symbolics.Unknown()
+            for i in eachindex(k)
+                new_varmap[k[i]] = v[i]
+                new_varmap[toterm(k[i])] = v[i]
+            end
+        end
+    end
+    return new_varmap
+end

ext/CatalystStructuralIdentifiabilityExtension/structural_identifiability_extension.jl

@@ -169,7 +169,7 @@ function make_osys(rs::ReactionSystem; remove_conserved = true)
         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))
+    osys = complete(make_rre_ode(rs; remove_conserved))
     vars = [unknowns(rs); parameters(rs)]
 
     # Computes equations for system conservation laws.

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, merge, check_units,
                         get_unit, check_equations, iscomplete
 
 import Base: (==), hash, size, getindex, setindex, isless, Sort.defalg, length, show
@@ -94,13 +94,14 @@ export isautonomous
 export reactionrates
 export isequivalent
 export set_default_noise_scaling
+export make_rre_ode, make_cle_sde, make_sck_jump, make_rre_algeqs
 
 # depreciated functions to remove in future releases
 export params, numparams
 
 # Conversions of the `ReactionSystem` structure.
 include("reactionsystem_conversions.jl")
-export ODEProblem, SDEProblem, JumpProblem, NonlinearProblem, DiscreteProblem,
+export ODEProblem, SDEProblem, JumpProblem, NonlinearProblem,
        SteadyStateProblem, JumpInputs
 export ismassaction, oderatelaw, jumpratelaw
 export symmap_to_varmap

src/dsl.jl

@@ -262,7 +262,6 @@ end
 
 # Function for creating a ReactionSystem structure (used by the @reaction_network macro).
 function make_reaction_system(ex::Expr, name)
-
     # Handle interpolation of variables in the input.
     ex = esc_dollars!(ex)
 
@@ -988,9 +987,8 @@ function recursive_escape_functions!(expr::ExprValues, syms_skip = [])
     (typeof(expr) != Expr) && (return expr)
     foreach(i -> expr.args[i] = recursive_escape_functions!(expr.args[i], syms_skip),
         1:length(expr.args))
-    if (expr.head == :call) && (expr.args[1] isa Symbol) &&
-       !isdefined(Catalyst, expr.args[1]) &&
-       expr.args[1] ∉ syms_skip
+    if (expr.head == :call) && (expr.args[1] isa Symbol) &&!isdefined(Catalyst, expr.args[1]) &&
+            expr.args[1] ∉ syms_skip
         expr.args[1] = esc(expr.args[1])
     end
     expr

src/latexify_recipes.jl

@@ -22,7 +22,7 @@ const LATEX_DEFS = CatalystLatexParams()
         return convert(ODESystem, rs)
     elseif form == :sde      # Returns SDE system code.
         mult_symbol --> ""
-        return convert(SDESystem, rs)
+        return make_cle_sde(rs)
     end
     error("Unrecognised form argument given: $form. This should be either reactions (default), :ode, or :sde.")
 end

src/reactionsystem.jl

@@ -275,8 +275,7 @@ Notes:
   the same units, and all reactions have rate laws with units of (species units) / (time
   units). Unit checking can be disabled by passing the keyword argument `checks=false`.
 """
-struct ReactionSystem{V <: NetworkProperties} <:
-       MT.AbstractTimeDependentSystem
+struct ReactionSystem{V <: NetworkProperties} <: MT.AbstractSystem
     """The equations (reactions and algebraic/differential) defining the system."""
     eqs::Vector{CatalystEqType}
     """The Reactions defining the system. """
@@ -373,7 +372,6 @@ struct ReactionSystem{V <: NetworkProperties} <:
             (hasnode(is_species_diff, eq.lhs) || hasnode(is_species_diff, eq.rhs)) &&
                 error("An equation ($eq) contains a differential with respect to a species. This is currently not supported. If this is a functionality you require, please raise an issue on the Catalyst GitHub page and we can consider the best way to implement it.")
         end
-
         rs = new{typeof(nps)}(
             eqs, rxs, iv, sivs, unknowns, spcs, ps, var_to_name, observed,
             name, systems, defaults, connection_type, nps, cls, cevs,
@@ -398,7 +396,7 @@ function ReactionSystem(eqs, iv, unknowns, ps;
         name = nothing,
         default_u0 = Dict(),
         default_p = Dict(),
-        defaults = _merge(Dict(default_u0), Dict(default_p)),
+        defaults = merge(Dict(default_u0), Dict(default_p)),
         connection_type = nothing,
         checks = true,
         networkproperties = nothing,
@@ -472,22 +470,24 @@ function ReactionSystem(eqs, iv, unknowns, ps;
     MT.process_variables!(var_to_name, defaults, unknowns′)
     MT.process_variables!(var_to_name, defaults, ps′)
     MT.collect_var_to_name!(var_to_name, eq.lhs for eq in observed)
-    #
+
     # Computes network properties.
     nps = if networkproperties === nothing
         NetworkProperties{Int, get_speciestype(iv′, unknowns′, systems)}()
     else
         networkproperties
     end
 
-    # Creates the continuous and discrete callbacks.
-    ccallbacks = MT.SymbolicContinuousCallbacks(continuous_events)
-    dcallbacks = MT.SymbolicDiscreteCallbacks(discrete_events)
+    # Creates the continuous and discrete events.
+    continuous_events, discrete_events = MT.create_symbolic_events(continuous_events, discrete_events)
+
+    # handles system metadata.
+    metadata === nothing ? Base.ImmutableDict{Symbol,Any}() : metadata
 
     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).
@@ -1385,20 +1385,6 @@ function make_empty_network(; iv = DEFAULT_IV, name = gensym(:ReactionSystem))
     ReactionSystem(Reaction[], iv, [], []; name = name)
 end
 
-# A helper function used in `flatten`.
-function getsubsystypes!(typeset::Set{Type}, sys::T) where {T <: MT.AbstractSystem}
-    push!(typeset, T)
-    for subsys in get_systems(sys)
-        getsubsystypes!(typeset, subsys)
-    end
-    typeset
-end
-
-function getsubsystypes(sys)
-    typeset = Set{Type}()
-    getsubsystypes!(typeset, sys)
-    typeset
-end
 
 """
     ModelingToolkit.flatten(rs::ReactionSystem)
@@ -1419,11 +1405,10 @@ Notes:
 function MT.flatten(rs::ReactionSystem; name = nameof(rs))
     isempty(get_systems(rs)) && return rs
 
-    # right now only NonlinearSystems and ODESystems can be handled as subsystems
-    subsys_types = getsubsystypes(rs)
+    # right now we only guarantee tht certain types of systems work with flatten
     allowed_types = (ReactionSystem, NonlinearSystem, ODESystem)
-    all(T -> any(T .<: allowed_types), subsys_types) ||
-        error("flattening is currently only supported for subsystems mixing ReactionSystems, NonlinearSystems and ODESystems.")
+    isnothing(get_systems(rs)) || all(is_allowed_subsystem, get_systems(rs)) ||
+        error("flattening is currently only supported for subsystems mixing ReactionSystems, and Systems withour noise equations and jumps.")
 
     ReactionSystem(equations(rs), get_iv(rs), unknowns(rs), parameters(rs);
         observed = MT.observed(rs),
@@ -1438,6 +1423,15 @@ function MT.flatten(rs::ReactionSystem; name = nameof(rs))
         metadata = MT.get_metadata(rs))
 end
 
+# Checks if a system is an allowed subsystem (i.e. no SDE parts and no jump).
+is_allowed_subsystem(sys::ReactionSystem) = true
+function is_allowed_subsystem(sys::System)
+    return (isnothing(MT.get_noise_eqs(sys)) || isempty(MT.get_noise_eqs(sys))) &&
+        (isnothing(MT.get_jumps(sys)) || isempty(MT.get_jumps(sys)))
+end
+# If neither a `ReactionSystem` or a `System`, it is something weird we do not know what it is.
+is_allowed_subsystem(sys::MT.AbstractSystem) = false
+
 function complete_check(sys, method)
     if MT.iscomplete(sys)
         error("$method with one or more `ReactionSystem`s requires systems to not be marked complete, but system: $(MT.get_name(sys)) is marked complete.")
@@ -1611,3 +1605,9 @@ unitless_exp(u) = iscall(u) && (operation(u) == ^) && (arguments(u)[1] == 1)
 function unitless_symvar(sym)
     return (sym isa Symbolics.CallWithMetadata) || (ModelingToolkit.get_unit(sym) == 1)
 end
+
+
+### Unsorted ###
+
+# Function previously used by ModelingToolkit.
+MT.refreshed_metadata(::Nothing) = MT.MetadataT() # FIXME: Type piracy