SciML · ChrisRackauckas · May 26, 2025 · May 21, 2025 · May 21, 2025 · May 21, 2025
diff --git a/ext/MTKBifurcationKitExt.jl b/ext/MTKBifurcationKitExt.jl
@@ -113,11 +113,12 @@ function BifurcationKit.BifurcationProblem(nsys::System,
     J = jac ? ofun.jac : nothing
 
     # Converts the input state guess.
-    u0_bif_vals = ModelingToolkit.varmap_to_vars(u0_bif,
-        unknowns(nsys);
-        defaults = ModelingToolkit.get_defaults(nsys))
-    p_vals = ModelingToolkit.varmap_to_vars(
-        ps, parameters(nsys); defaults = ModelingToolkit.get_defaults(nsys))
+    u0_bif = ModelingToolkit.to_varmap(u0_bif, unknowns(nsys))
+    u0_buf = merge(ModelingToolkit.get_defaults(nsys), u0_bif)
+    u0_bif_vals = ModelingToolkit.varmap_to_vars(u0_bif, unknowns(nsys))
+    ps = ModelingToolkit.to_varmap(ps, parameters(nsys))
+    ps = merge(ModelingToolkit.get_defaults(nsys), ps)
+    p_vals = ModelingToolkit.varmap_to_vars(ps, parameters(nsys))
 
     # Computes bifurcation parameter and the plotting function.
     bif_idx = findfirst(isequal(bif_par), parameters(nsys))

diff --git a/src/ModelingToolkit.jl b/src/ModelingToolkit.jl
@@ -128,10 +128,6 @@ $(TYPEDEF)
 TODO
 """
 abstract type AbstractSystem end
-abstract type AbstractTimeDependentSystem <: AbstractSystem end
-abstract type AbstractTimeIndependentSystem <: AbstractSystem end
-abstract type AbstractMultivariateSystem <: AbstractSystem end
-abstract type AbstractOptimizationSystem <: AbstractTimeIndependentSystem end
 
 function independent_variable end
 
@@ -150,7 +146,6 @@ include("independent_variables.jl")
 include("constants.jl")
 
 include("utils.jl")
-include("domains.jl")
 
 include("systems/index_cache.jl")
 include("systems/parameter_buffer.jl")
@@ -265,25 +260,19 @@ PrecompileTools.@compile_workload begin
     end
 end
 
-export AbstractTimeDependentSystem,
-       AbstractTimeIndependentSystem,
-       AbstractMultivariateSystem
-
-export ODEFunction, ODEFunctionExpr, ODEProblemExpr, convert_system_indepvar,
+export ODEFunction, convert_system_indepvar,
        System, OptimizationSystem, JumpSystem, SDESystem, NonlinearSystem
-export DAEFunctionExpr, DAEProblemExpr
-export SDEFunction, SDEFunctionExpr, SDEProblemExpr
+export SDEFunction
 export SystemStructure
-export DiscreteProblem, DiscreteFunction, DiscreteFunctionExpr
-export ImplicitDiscreteProblem, ImplicitDiscreteFunction,
-       ImplicitDiscreteFunctionExpr
+export DiscreteProblem, DiscreteFunction
+export ImplicitDiscreteProblem, ImplicitDiscreteFunction
 export ODEProblem, SDEProblem
-export NonlinearFunction, NonlinearFunctionExpr
-export NonlinearProblem, NonlinearProblemExpr
-export IntervalNonlinearFunction, IntervalNonlinearFunctionExpr
-export IntervalNonlinearProblem, IntervalNonlinearProblemExpr
-export OptimizationProblem, OptimizationProblemExpr, constraints
-export SteadyStateProblem, SteadyStateProblemExpr
+export NonlinearFunction
+export NonlinearProblem
+export IntervalNonlinearFunction
+export IntervalNonlinearProblem
+export OptimizationProblem, constraints
+export SteadyStateProblem
 export JumpProblem
 export alias_elimination, flatten
 export connect, domain_connect, @connector, Connection, AnalysisPoint, Flow, Stream,
@@ -315,7 +304,7 @@ export calculate_jacobian, generate_jacobian, generate_rhs, generate_custom_func
 export calculate_control_jacobian, generate_control_jacobian
 export calculate_tgrad, generate_tgrad
 export generate_cost, calculate_cost_gradient, generate_cost_gradient
-export calculate_factorized_W, generate_factorized_W
+export calculate_factorized_W
 export calculate_cost_hessian, generate_cost_hessian
 export calculate_massmatrix, generate_diffusion_function
 export stochastic_integral_transform

diff --git a/src/domains.jl b/src/domains.jl
diff --git a/src/problems/optimizationproblem.jl b/src/problems/optimizationproblem.jl
@@ -114,9 +114,14 @@ function SciMLBase.OptimizationProblem{iip}(
     end
 
     ps = parameters(sys)
-    defs = merge(defaults(sys), to_varmap(op, dvs))
-    lb = varmap_to_vars(dvs .=> lb, dvs; defaults = defs, tofloat = false)
-    ub = varmap_to_vars(dvs .=> ub, dvs; defaults = defs, tofloat = false)
+    defs = defaults(sys)
+    op = to_varmap(op, dvs)
+    lbmap = merge(op, AnyDict(dvs .=> lb))
+    _, _ = build_operating_point!(sys, lbmap, Dict(), Dict(), defs, dvs, ps)
+    lb = varmap_to_vars(lbmap, dvs; tofloat = false)
+    ubmap = merge(op, AnyDict(dvs .=> ub))
+    _, _ = build_operating_point!(sys, ubmap, Dict(), Dict(), defs, dvs, ps)
+    ub = varmap_to_vars(ubmap, dvs; tofloat = false)
 
     if !isnothing(lb) && all(lb .== -Inf) && !isnothing(ub) && all(ub .== Inf)
         lb = nothing

diff --git a/src/structural_transformation/StructuralTransformations.jl b/src/structural_transformation/StructuralTransformations.jl
@@ -21,9 +21,8 @@ using ModelingToolkit: System, AbstractSystem, var_from_nested_derivative, Diffe
                        has_tearing_state, defaults, InvalidSystemException,
                        ExtraEquationsSystemException,
                        ExtraVariablesSystemException,
-                       vars!,
+                       vars!, invalidate_cache!,
                        IncrementalCycleTracker, add_edge_checked!, topological_sort,
-                       invalidate_cache!, Substitutions, get_or_construct_tearing_state,
                        filter_kwargs, lower_varname_with_unit,
                        lower_shift_varname_with_unit, setio, SparseMatrixCLIL,
                        get_fullvars, has_equations, observed,

diff --git a/src/structural_transformation/symbolics_tearing.jl b/src/structural_transformation/symbolics_tearing.jl
@@ -734,7 +734,6 @@ function update_simplified_system!(
 
     @set! sys.eqs = neweqs
     @set! sys.observed = obs
-    # @set! sys.substitutions = Substitutions(subeqs, deps)
 
     # Only makes sense for time-dependent
     if ModelingToolkit.has_schedule(sys)

diff --git a/src/structural_transformation/utils.jl b/src/structural_transformation/utils.jl
@@ -352,8 +352,14 @@ function but_ordered_incidence(ts::TearingState, varmask = highest_order_variabl
     mm[[var_eq_matching[v] for v in vordering if var_eq_matching[v] isa Int], vordering], bb
 end
 
-# debugging use
-function reordered_matrix(sys, torn_matching)
+"""
+    $(TYPEDSIGNATURES)
+
+Given a system `sys` and torn variable-equation matching `torn_matching`, return the sparse
+incidence matrix of the system with SCCs grouped together, and each SCC sorted to contain
+the analytically solved equations/variables before the unsolved ones.
+"""
+function reordered_matrix(sys::System, torn_matching)
     s = TearingState(sys)
     complete!(s.structure)
     @unpack graph = s.structure
@@ -432,23 +438,6 @@ function torn_system_jacobian_sparsity(sys)
     return sparse(I, J, true, neqs, neqs)
 end
 
-###
-### Nonlinear equation(s) solving
-###
-
-@noinline function nlsolve_failure(rc)
-    error("The nonlinear solver failed with the return code $rc.")
-end
-
-function numerical_nlsolve(f, u0, p)
-    prob = NonlinearProblem{false}(f, u0, p)
-    sol = solve(prob, SimpleNewtonRaphson())
-    rc = sol.retcode
-    (rc == ReturnCode.Success) || nlsolve_failure(rc)
-    # TODO: robust initial guess, better debugging info, and residual check
-    sol.u
-end
-
 ###
 ### Misc
 ###