[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:

.gitignore

@@ -10,4 +10,8 @@ Manifest.toml
 .vscode
 .vscode/*
 
-LocalPreferences.toml
+LocalPreferences.toml
+
+# claude
+.claude
+.claude/*

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/Project.toml

@@ -71,7 +71,7 @@ ModelingToolkit = "9.69"
 NetworkLayout = "0.4"
 NonlinearSolve = "3.12, 4"
 Optim = "1.9"
-Optimization = "4"
+Optimization = "4, 5"
 OptimizationBBO = "0.4"
 OptimizationEvolutionary = "0.4"
 OptimizationNLopt = "0.3"

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