Merge branch 'master' into BVP-with-constraints

Author: vyudu

.github/workflows/Downstream.yml

@@ -37,6 +37,7 @@ jobs:
           - {user: SciML, repo: MethodOfLines.jl, group: Interface}
           - {user: SciML, repo: MethodOfLines.jl, group: 2D_Diffusion}
           - {user: SciML, repo: MethodOfLines.jl, group: DAE}
+          - {user: SciML, repo: ModelingToolkitNeuralNets.jl, group: All}
     steps:
       - uses: actions/checkout@v4
       - uses: julia-actions/setup-julia@v1

Project.toml

@@ -121,6 +121,7 @@ LabelledArrays = "1.3"
 Latexify = "0.11, 0.12, 0.13, 0.14, 0.15, 0.16"
 Libdl = "1"
 LinearAlgebra = "1"
+Logging = "1"
 MLStyle = "0.4.17"
 ModelingToolkitStandardLibrary = "2.19"
 NaNMath = "0.3, 1"
@@ -130,7 +131,7 @@ OrderedCollections = "1"
 OrdinaryDiffEq = "6.82.0"
 OrdinaryDiffEqCore = "1.15.0"
 OrdinaryDiffEqDefault = "1.2"
-OrdinaryDiffEqNonlinearSolve = "1.3.0"
+OrdinaryDiffEqNonlinearSolve = "1.5.0"
 PrecompileTools = "1"
 REPL = "1"
 RecursiveArrayTools = "3.26"
@@ -146,9 +147,10 @@ SparseArrays = "1"
 SpecialFunctions = "0.7, 0.8, 0.9, 0.10, 1.0, 2"
 StaticArrays = "0.10, 0.11, 0.12, 1.0"
 StochasticDelayDiffEq = "1.8.1"
+StochasticDiffEq = "6.72.1"
 SymbolicIndexingInterface = "0.3.37"
 SymbolicUtils = "3.14"
-Symbolics = "6.29"
+Symbolics = "6.29.1"
 URIs = "1"
 UnPack = "0.1, 1.0"
 Unitful = "1.1"
@@ -167,6 +169,7 @@ ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 Ipopt = "b6b21f68-93f8-5de0-b562-5493be1d77c9"
 Ipopt_jll = "9cc047cb-c261-5740-88fc-0cf96f7bdcc7"
 JET = "c3a54625-cd67-489e-a8e7-0a5a0ff4e31b"
+Logging = "56ddb016-857b-54e1-b83d-db4d58db5568"
 ModelingToolkitStandardLibrary = "16a59e39-deab-5bd0-87e4-056b12336739"
 NonlinearSolve = "8913a72c-1f9b-4ce2-8d82-65094dcecaec"
 Optimization = "7f7a1694-90dd-40f0-9382-eb1efda571ba"
@@ -190,4 +193,4 @@ Sundials = "c3572dad-4567-51f8-b174-8c6c989267f4"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["AmplNLWriter", "BenchmarkTools", "BoundaryValueDiffEqMIRK", "BoundaryValueDiffEqAscher", "ControlSystemsBase", "DataInterpolations", "DelayDiffEq", "NonlinearSolve", "ForwardDiff", "Ipopt", "Ipopt_jll", "ModelingToolkitStandardLibrary", "Optimization", "OptimizationOptimJL", "OptimizationMOI", "OrdinaryDiffEq", "OrdinaryDiffEqCore", "OrdinaryDiffEqDefault", "REPL", "Random", "ReferenceTests", "SafeTestsets", "StableRNGs", "Statistics", "SteadyStateDiffEq", "Test", "StochasticDiffEq", "Sundials", "StochasticDelayDiffEq", "Pkg", "JET", "OrdinaryDiffEqNonlinearSolve"]
+test = ["AmplNLWriter", "BenchmarkTools", "BoundaryValueDiffEqMIRK", "BoundaryValueDiffEqAscher", "ControlSystemsBase", "DataInterpolations", "DelayDiffEq", "NonlinearSolve", "ForwardDiff", "Ipopt", "Ipopt_jll", "ModelingToolkitStandardLibrary", "Optimization", "OptimizationOptimJL", "OptimizationMOI", "OrdinaryDiffEq", "OrdinaryDiffEqCore", "OrdinaryDiffEqDefault", "REPL", "Random", "ReferenceTests", "SafeTestsets", "StableRNGs", "Statistics", "SteadyStateDiffEq", "Test", "StochasticDiffEq", "Sundials", "StochasticDelayDiffEq", "Pkg", "JET", "OrdinaryDiffEqNonlinearSolve", "Logging"]

ext/MTKBifurcationKitExt.jl

@@ -97,7 +97,10 @@ function BifurcationKit.BifurcationProblem(nsys::NonlinearSystem,
     @set! nsys.index_cache = nothing # force usage of a parameter vector instead of `MTKParameters`
     # Creates F and J functions.
     ofun = NonlinearFunction(nsys; jac = jac)
-    F = ofun.f
+    F = let f = ofun.f
+        _f(resid, u, p) = (f(resid, u, p); resid)
+        _f(u, p) = f(u, p)
+    end
     J = jac ? ofun.jac : nothing
 
     # Converts the input state guess.
@@ -136,6 +139,7 @@ function BifurcationKit.BifurcationProblem(nsys::NonlinearSystem,
         args...;
         record_from_solution = record_from_solution,
         J = J,
+        inplace = true,
         kwargs...)
 end
 

ext/MTKChainRulesCoreExt.jl

@@ -16,7 +16,12 @@ end
 
 function subset_idxs(idxs, portion, template)
     ntuple(Val(length(template))) do subi
-        [Base.tail(idx.idx) for idx in idxs if idx.portion == portion && idx.idx[1] == subi]
+        result = [Base.tail(idx.idx)
+                  for idx in idxs if idx.portion == portion && idx.idx[1] == subi]
+        if isempty(result)
+            result = []
+        end
+        result
     end
 end
 

src/linearization.jl

@@ -170,8 +170,8 @@ function (linfun::LinearizationFunction)(u, p, t)
     if u !== nothing # Handle systems without unknowns
         linfun.num_states == length(u) ||
             error("Number of unknown variables ($(linfun.num_states)) does not match the number of input unknowns ($(length(u)))")
-        integ_cache = linfun.caches
-        integ = MockIntegrator{true}(u, p, t, integ_cache)
+        integ_cache = (linfun.caches,)
+        integ = MockIntegrator{true}(u, p, t, integ_cache, nothing)
         u, p, success = SciMLBase.get_initial_values(
             linfun.prob, integ, fun, linfun.initializealg, Val(true);
             linfun.initialize_kwargs...)
@@ -218,7 +218,7 @@ Mock `DEIntegrator` to allow using `CheckInit` without having to create a new in
 
 $(TYPEDFIELDS)
 """
-struct MockIntegrator{iip, U, P, T, C} <: SciMLBase.DEIntegrator{Nothing, iip, U, T}
+struct MockIntegrator{iip, U, P, T, C, O} <: SciMLBase.DEIntegrator{Nothing, iip, U, T}
     """
     The state vector.
     """
@@ -235,10 +235,14 @@ struct MockIntegrator{iip, U, P, T, C} <: SciMLBase.DEIntegrator{Nothing, iip, U
     The integrator cache.
     """
     cache::C
+    """
+    Integrator "options" for `CheckInit`.
+    """
+    opts::O
 end
 
-function MockIntegrator{iip}(u::U, p::P, t::T, cache::C) where {iip, U, P, T, C}
-    return MockIntegrator{iip, U, P, T, C}(u, p, t, cache)
+function MockIntegrator{iip}(u::U, p::P, t::T, cache::C, opts::O) where {iip, U, P, T, C, O}
+    return MockIntegrator{iip, U, P, T, C, O}(u, p, t, cache, opts)
 end
 
 SymbolicIndexingInterface.state_values(integ::MockIntegrator) = integ.u

src/structural_transformation/symbolics_tearing.jl

@@ -647,6 +647,7 @@ function cse_and_array_hacks(sys, obs, subeqs, unknowns, neweqs; cse = true, arr
         # HACK 1
         if cse && is_getindexed_array(rhs)
             rhs_arr = arguments(rhs)[1]
+            iscall(rhs_arr) && operation(rhs_arr) isa Symbolics.Operator && continue
             if !haskey(rhs_to_tempvar, rhs_arr)
                 tempvar = gensym(Symbol(lhs))
                 N = length(rhs_arr)
@@ -719,6 +720,7 @@ function cse_and_array_hacks(sys, obs, subeqs, unknowns, neweqs; cse = true, arr
         Symbolics.shape(sym) != Symbolics.Unknown() || continue
         arg1 = arguments(sym)[1]
         cnt = get(arr_obs_occurrences, arg1, 0)
+        cnt == 0 && continue
         arr_obs_occurrences[arg1] = cnt + 1
     end
     for eq in neweqs

src/systems/abstractsystem.jl

@@ -669,7 +669,11 @@ end
 
 # This is required so `fast_substitute` works
 function SymbolicUtils.maketerm(::Type{<:BasicSymbolic}, ::Initial, args, meta)
-    return metadata(Initial()(args...), meta)
+    val = Initial()(args...)
+    if symbolic_type(val) == NotSymbolic()
+        return val
+    end
+    return metadata(val, meta)
 end
 
 function add_initialization_parameters(sys::AbstractSystem)

src/systems/codegen_utils.jl

@@ -12,8 +12,14 @@ end
 
 Given the arguments to `build_function_wrapper`, return a list of assignments which
 reconstruct array variables if they are present scalarized in `args`.
+
+# Keyword Arguments
+
+- `argument_name` a function of the form `(::Int) -> Symbol` which takes the index of
+  an argument to the generated function and returns the name of the argument in the
+  generated function.
 """
-function array_variable_assignments(args...)
+function array_variable_assignments(args...; argument_name = generated_argument_name)
     # map array symbolic to an identically sized array where each element is (buffer_idx, idx_in_buffer)
     var_to_arridxs = Dict{BasicSymbolic, Array{Tuple{Int, Int}}}()
     for (i, arg) in enumerate(args)
@@ -60,12 +66,12 @@ function array_variable_assignments(args...)
             end
             # view and reshape
 
-            expr = term(reshape, term(view, generated_argument_name(buffer_idx), idxs),
+            expr = term(reshape, term(view, argument_name(buffer_idx), idxs),
                 size(arrvar))
         else
             elems = map(idxs) do idx
                 i, j = idx
-                term(getindex, generated_argument_name(i), j)
+                term(getindex, argument_name(i), j)
             end
             # use `MakeArray` syntax and generate a stack-allocated array
             expr = term(SymbolicUtils.Code.create_array, SArray, nothing,

src/systems/diffeqs/abstractodesystem.jl

@@ -1455,7 +1455,9 @@ function InitializationProblem{iip, specialize}(sys::AbstractSystem,
         guesses = Dict()
     end
 
-    u0map = merge(ModelingToolkit.guesses(sys), todict(guesses), todict(u0map))
+    filter_missing_values!(u0map)
+    filter_missing_values!(parammap)
+    u0map = merge(ModelingToolkit.guesses(sys), todict(guesses), u0map)
 
     fullmap = merge(u0map, parammap)
     u0T = Union{}

src/systems/diffeqs/odesystem.jl

@@ -505,11 +505,16 @@ function build_explicit_observed_function(sys, ts;
         end
     end
     allsyms = Set(all_symbols(sys))
+    iv = has_iv(sys) ? get_iv(sys) : nothing
     for var in vs
         var = unwrap(var)
         newvar = get(ns_map, var, nothing)
         if newvar !== nothing
             namespace_subs[var] = newvar
+            var = newvar
+        end
+        if throw && !var_in_varlist(var, allsyms, iv)
+            Base.throw(ArgumentError("Symbol $var is not present in the system."))
         end
     end
     ts = fast_substitute(ts, namespace_subs)
@@ -551,12 +556,14 @@ function build_explicit_observed_function(sys, ts;
     if fns isa Tuple
         oop, iip = eval_or_rgf.(fns; eval_expression, eval_module)
         f = GeneratedFunctionWrapper{(
-            p_start, length(args) - length(ps) + 1, is_split(sys))}(oop, iip)
+            p_start + is_dde(sys), length(args) - length(ps) + 1 + is_dde(sys), is_split(sys))}(
+            oop, iip)
         return return_inplace ? (f, f) : f
     else
         f = eval_or_rgf(fns; eval_expression, eval_module)
         f = GeneratedFunctionWrapper{(
-            p_start, length(args) - length(ps) + 1, is_split(sys))}(f, nothing)
+            p_start + is_dde(sys), length(args) - length(ps) + 1 + is_dde(sys), is_split(sys))}(
+            f, nothing)
         return f
     end
 end