Merge branch 'SciML:master' into func-affect

Author: dodoplus

Project.toml

@@ -1,7 +1,7 @@
 name = "ModelingToolkit"
 uuid = "961ee093-0014-501f-94e3-6117800e7a78"
 authors = ["Chris Rackauckas <[email protected]>"]
-version = "8.14.1"
+version = "8.15.1"
 
 [deps]
 AbstractTrees = "1520ce14-60c1-5f80-bbc7-55ef81b5835c"
@@ -11,7 +11,6 @@ ConstructionBase = "187b0558-2788-49d3-abe0-74a17ed4e7c9"
 DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
 DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e"
 DiffEqCallbacks = "459566f4-90b8-5000-8ac3-15dfb0a30def"
-DiffEqJump = "c894b116-72e5-5b58-be3c-e6d8d4ac2b12"
 DiffRules = "b552c78f-8df3-52c6-915a-8e097449b14b"
 Distributed = "8ba89e20-285c-5b6f-9357-94700520ee1b"
 Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
@@ -21,6 +20,7 @@ Graphs = "86223c79-3864-5bf0-83f7-82e725a168b6"
 IfElse = "615f187c-cbe4-4ef1-ba3b-2fcf58d6d173"
 InteractiveUtils = "b77e0a4c-d291-57a0-90e8-8db25a27a240"
 JuliaFormatter = "98e50ef6-434e-11e9-1051-2b60c6c9e899"
+JumpProcesses = "ccbc3e58-028d-4f4c-8cd5-9ae44345cda5"
 LabelledArrays = "2ee39098-c373-598a-b85f-a56591580800"
 Latexify = "23fbe1c1-3f47-55db-b15f-69d7ec21a316"
 Libdl = "8f399da3-3557-5675-b5ff-fb832c97cbdb"
@@ -50,14 +50,14 @@ ConstructionBase = "1"
 DataStructures = "0.17, 0.18"
 DiffEqBase = "6.83.0"
 DiffEqCallbacks = "2.16"
-DiffEqJump = "7.0, 8"
 DiffRules = "0.1, 1.0"
 Distributions = "0.23, 0.24, 0.25"
 DocStringExtensions = "0.7, 0.8, 0.9"
 DomainSets = "0.5"
 Graphs = "1.5.2"
 IfElse = "0.1"
 JuliaFormatter = "1"
+JumpProcesses = "9"
 LabelledArrays = "1.3"
 Latexify = "0.11, 0.12, 0.13, 0.14, 0.15"
 MacroTools = "0.5"

src/ModelingToolkit.jl

@@ -12,7 +12,7 @@ using Latexify, Unitful, ArrayInterfaceCore
 using MacroTools
 @reexport using UnPack
 using Setfield, ConstructionBase
-using DiffEqJump
+using JumpProcesses
 using DataStructures
 using SpecialFunctions, NaNMath
 using RuntimeGeneratedFunctions

src/structural_transformation/bipartite_tearing/modia_tearing.jl

@@ -37,57 +37,17 @@ end
 
 function tear_graph_modia(structure::SystemStructure; varfilter = v -> true,
                           eqfilter = eq -> true)
-    # It would be possible here to simply iterate over all variables and attempt to
-    # use tearEquations! to produce a matching that greedily selects the minimal
-    # number of torn variables. However, we can do this process faster if we first
-    # compute the strongly connected components. In the absence of cycles and
-    # non-solvability, a maximal matching on the original graph will give us an
-    # optimal assignment. However, even with cycles, we can use the maximal matching
-    # to give us a good starting point for a good matching and then proceed to
-    # reverse edges in each scc to improve the solution. Note that it is possible
-    # to have optimal solutions that cannot be found by this process. We will not
-    # find them here [TODO: It would be good to have an explicit example of this.]
-
     @unpack graph, solvable_graph = structure
     var_eq_matching = complete(maximal_matching(graph, eqfilter, varfilter))
     var_sccs::Vector{Union{Vector{Int}, Int}} = find_var_sccs(graph, var_eq_matching)
 
-    # Here, we're using a maximal matching on the post-pantelides system to find
-    # the strongly connected components of the system (of variables that depend
-    # on each other). The strongly connected components are unique, however, the
-    # maximal matching itself is not. Every maximal matching gives rise to the
-    # same set of strongly connected components, but the associated equations need
-    # not be the same. In the absence of solvability constraints, this may be a
-    # small issue, but here it is possible that an equation got assigned to an
-    # scc that cannot actually use it for solving a variable, but still precludes
-    # another scc from using it. To avoid this, we delete any assignments that
-    # are not in the solvable graph and extend the set of considered eqauations
-    # below.
-    for var in ndsts(solvable_graph)
-        var_eq_matching[var] === unassigned && continue
-        if !(BipartiteEdge(var, var_eq_matching[var]) in solvable_graph)
-            var_eq_matching[var] = unassigned
-        end
-    end
-
     for vars in var_sccs
         filtered_vars = filter(varfilter, vars)
         ieqs = Int[var_eq_matching[v]
                    for v in filtered_vars if var_eq_matching[v] !== unassigned]
         for var in vars
             var_eq_matching[var] = unassigned
         end
-        for var in filtered_vars
-            # Add any equations that we may not have been able to use earlier to see
-            # if a different matching may have been possible.
-            for eq′ in 𝑑neighbors(solvable_graph, var)
-                eqfilter(eq′) || continue
-                eq′ in ieqs && continue
-                if invview(var_eq_matching)[eq′] === unassigned
-                    push!(ieqs, eq′)
-                end
-            end
-        end
         tear_graph_block_modia!(var_eq_matching, graph, solvable_graph, ieqs, filtered_vars)
     end
 

src/systems/alias_elimination.jl

@@ -57,13 +57,7 @@ function alias_elimination(sys)
 
     newstates = []
     for j in eachindex(fullvars)
-        if j in keys(ag)
-            # Put back equations for alias eliminated dervars
-            if isdervar(state.structure, j) &&
-               !(invview(state.structure.var_to_diff)[j] in keys(ag))
-                push!(eqs, fullvars[j] ~ subs[fullvars[j]])
-            end
-        else
+        if !(j in keys(ag))
             isdervar(state.structure, j) || push!(newstates, fullvars[j])
         end
     end
@@ -198,7 +192,6 @@ struct AliasGraphKeySet <: AbstractSet{Int}
 end
 Base.keys(ag::AliasGraph) = AliasGraphKeySet(ag)
 Base.iterate(agk::AliasGraphKeySet, state...) = Base.iterate(agk.ag.eliminated, state...)
-Base.length(agk::AliasGraphKeySet) = Base.length(agk.ag.eliminated)
 function Base.in(i::Int, agk::AliasGraphKeySet)
     aliasto = agk.ag.aliasto
     1 <= i <= length(aliasto) && aliasto[i] !== nothing
@@ -217,11 +210,9 @@ function aag_bareiss!(graph, var_to_diff, mm_orig::SparseMatrixCLIL)
         is_linear_equations[e] = true
     end
 
-    # For now, only consider variables linear that are not differentiated.
-    # We could potentially apply the same logic to variables whose derivative
-    # is also linear, but that's a TODO.
-    diff_to_var = invview(var_to_diff)
-    is_linear_variables = .&(isnothing.(var_to_diff), isnothing.(diff_to_var))
+    # Variables that are highest order differentiated cannot be states of an ODE
+    is_not_potential_state = isnothing.(var_to_diff)
+    is_linear_variables = copy(is_not_potential_state)
     for i in 𝑠vertices(graph)
         is_linear_equations[i] && continue
         for j in 𝑠neighbors(graph, i)
@@ -239,9 +230,11 @@ function aag_bareiss!(graph, var_to_diff, mm_orig::SparseMatrixCLIL)
                 r !== nothing && return r
                 rank1 = k - 1
             end
-            # TODO: It would be better to sort the variables by
-            # derivative order here to enable more elimination
-            # opportunities.
+            if rank2 === nothing
+                r = find_masked_pivot(is_not_potential_state, M, k)
+                r !== nothing && return r
+                rank2 = k - 1
+            end
             return find_masked_pivot(nothing, M, k)
         end
         function find_and_record_pivot(M, k)
@@ -256,9 +249,10 @@ function aag_bareiss!(graph, var_to_diff, mm_orig::SparseMatrixCLIL)
         end
         bareiss_ops = ((M, i, j) -> nothing, myswaprows!,
                        bareiss_update_virtual_colswap_mtk!, bareiss_zero!)
-        rank2, = bareiss!(M, bareiss_ops; find_pivot = find_and_record_pivot)
-        rank1 = something(rank1, rank2)
-        (rank1, rank2, pivots)
+        rank3, = bareiss!(M, bareiss_ops; find_pivot = find_and_record_pivot)
+        rank1 = something(rank1, rank3)
+        rank2 = something(rank2, rank3)
+        (rank1, rank2, rank3, pivots)
     end
 
     return mm, solvable_variables, do_bareiss!(mm, mm_orig)
@@ -272,27 +266,16 @@ function alias_eliminate_graph!(graph, var_to_diff, mm_orig::SparseMatrixCLIL)
     # variables`.
     #
     # `do_bareiss` conceptually gives us this system:
-    # rank1 | [ M₁₁  M₁₂ | M₁₃ ]   [v₁] = [0]
-    # rank2 | [ 0    M₂₂ | M₂₃ ] P [v₂] = [0]
+    # rank1 | [ M₁₁  M₁₂ | M₁₃  M₁₄ ]   [v₁] = [0]
+    # rank2 | [ 0    M₂₂ | M₂₃  M₂₄ ] P [v₂] = [0]
     # -------------------|------------------------
-    #         [ 0    0   | 0   ]   [v₃] = [0]
-    #
-    # Where `v₁` are the purely linear variables (i.e. those that only appear in linear equations),
-    # `v₂` are the variables that may be potentially solved by the linear system and v₃ are the variables
-    # that contribute to the equations, but are not solved by the linear system. Note
-    # that the complete system may be larger than the linear subsystem and include variables
-    # that do not appear here.
-    mm, solvable_variables, (rank1, rank2, pivots) = aag_bareiss!(graph, var_to_diff,
-                                                                  mm_orig)
+    # rank3 | [ 0    0   | M₃₃  M₃₄ ]   [v₃] = [0]
+    #         [ 0    0   | 0    0   ]   [v₄] = [0]
+    mm, solvable_variables, (rank1, rank2, rank3, pivots) = aag_bareiss!(graph, var_to_diff,
+                                                                         mm_orig)
 
     # Step 2: Simplify the system using the Bareiss factorization
-
     ag = AliasGraph(size(mm, 2))
-
-    # First, eliminate variables that only appear in linear equations and were removed
-    # completely from the coefficient matrix. These are technically singularities in
-    # the matrix, but assigning them to 0 is a feasible assignment and works well in
-    # practice.
     for v in setdiff(solvable_variables, @view pivots[1:rank1])
         ag[v] = 0
     end
@@ -304,7 +287,11 @@ function alias_eliminate_graph!(graph, var_to_diff, mm_orig::SparseMatrixCLIL)
     function lss!(ei::Integer)
         vi = pivots[ei]
         may_eliminate = true
-        locally_structure_simplify!((@view mm[ei, :]), vi, ag, var_to_diff)
+        for v in 𝑠neighbors(graph, mm.nzrows[ei])
+            # the differentiated variable cannot be eliminated
+            may_eliminate &= isnothing(diff_to_var[v]) && isnothing(var_to_diff[v])
+        end
+        locally_structure_simplify!((@view mm[ei, :]), vi, ag, may_eliminate)
     end
 
     # Step 2.1: Go backwards, collecting eliminated variables and substituting
@@ -346,7 +333,7 @@ function exactdiv(a::Integer, b)
     return d
 end
 
-function locally_structure_simplify!(adj_row, pivot_col, ag, var_to_diff)
+function locally_structure_simplify!(adj_row, pivot_col, ag, may_eliminate)
     pivot_val = adj_row[pivot_col]
     iszero(pivot_val) && return false
 
@@ -388,36 +375,21 @@ function locally_structure_simplify!(adj_row, pivot_col, ag, var_to_diff)
         end
     end
 
-    if nirreducible <= 1
+    if may_eliminate && nirreducible <= 1
         # There were only one or two terms left in the equation (including the
         # pivot variable). We can eliminate the pivot variable.
         #
         # Note that when `nirreducible <= 1`, `alias_candidate` is uniquely
         # determined.
         if alias_candidate !== 0
-            # Verify that the derivative depth of the variable is at least
-            # as deep as that of the alias, otherwise, we can't eliminate.
-            pivot_var = pivot_col
-            alias_var = alias_candidate[2]
-            while (pivot_var = var_to_diff[pivot_col]) !== nothing
-                alias_var = var_to_diff[alias_var]
-                alias_var === nothing && return false
-            end
             d, r = divrem(alias_candidate[1], pivot_val)
             if r == 0 && (d == 1 || d == -1)
                 alias_candidate = -d => alias_candidate[2]
             else
                 return false
             end
         end
-        diff_alias_candidate(ac) = ac === 0 ? 0 : ac[1] => var_to_diff[ac[2]]
-        while true
-            @assert !haskey(ag, pivot_col)
-            ag[pivot_col] = alias_candidate
-            pivot_col = var_to_diff[pivot_col]
-            pivot_col === nothing && break
-            alias_candidate = diff_alias_candidate(alias_candidate)
-        end
+        ag[pivot_col] = alias_candidate
         zero!(adj_row)
         return true
     end

src/systems/jumps/jumpsystem.jl

@@ -212,7 +212,7 @@ and no SDEs associated with the system.
 
 Continuing the example from the [`JumpSystem`](@ref) definition:
 ```julia
-using DiffEqBase, DiffEqJump
+using DiffEqBase, JumpProcesses
 u₀map = [S => 999, I => 1, R => 0]
 parammap = [β => .1/1000, γ => .01]
 tspan = (0.0, 250.0)
@@ -264,7 +264,7 @@ and no SDEs associated with the system.
 
 Continuing the example from the [`JumpSystem`](@ref) definition:
 ```julia
-using DiffEqBase, DiffEqJump
+using DiffEqBase, JumpProcesses
 u₀map = [S => 999, I => 1, R => 0]
 parammap = [β => .1/1000, γ => .01]
 tspan = (0.0, 250.0)
@@ -305,7 +305,7 @@ jprob = JumpProblem(js, dprob, Direct())
 sol = solve(jprob, SSAStepper())
 ```
 """
-function DiffEqJump.JumpProblem(js::JumpSystem, prob, aggregator; kwargs...)
+function JumpProcesses.JumpProblem(js::JumpSystem, prob, aggregator; kwargs...)
     statetoid = Dict(value(state) => i for (i, state) in enumerate(states(js)))
     eqs = equations(js)
     invttype = prob.tspan[1] === nothing ? Float64 : typeof(1 / prob.tspan[2])
@@ -421,6 +421,6 @@ function (ratemap::JumpSysMajParamMapper{U, V, W})(maj::MassActionJump, newparam
                                       value(substitute(ratemap.paramexprs[i],
                                                        ratemap.subdict)))
     end
-    scale_rates && DiffEqJump.scalerates!(maj.scaled_rates, maj.reactant_stoch)
+    scale_rates && JumpProcesses.scalerates!(maj.scaled_rates, maj.reactant_stoch)
     nothing
 end

src/systems/systemstructure.jl

@@ -348,6 +348,7 @@ function linear_subsys_adjmat(state::TransformationState)
     cadj = Vector{Int}[]
     coeffs = Int[]
     for (i, eq) in enumerate(eqs)
+        isdiffeq(eq) && continue
         empty!(coeffs)
         linear_term = 0
         all_int_vars = true

src/variables.jl

@@ -237,16 +237,17 @@ function tunable_parameters(sys, p = parameters(sys); default = false)
 end
 
 """
-    getbounds(sys::ModelingToolkit.AbstractSystem)
+    getbounds(sys::ModelingToolkit.AbstractSystem, p = parameters(sys))
 
 Returns a dict with pairs `p => (lb, ub)` mapping parameters of `sys` to lower and upper bounds.
 Create parameters with bounds like this
 ```
 @parameters p [bounds=(-1, 1)]
 ```
+
+To obtain state bounds, call `getbounds(sys, states(sys))`
 """
-function getbounds(sys::ModelingToolkit.AbstractSystem)
-    p = parameters(sys)
+function getbounds(sys::ModelingToolkit.AbstractSystem, p = parameters(sys))
     Dict(p .=> getbounds.(p))
 end
 

test/dep_graphs.jl

@@ -1,5 +1,5 @@
 using Test
-using ModelingToolkit, Graphs, DiffEqJump
+using ModelingToolkit, Graphs, JumpProcesses
 
 import ModelingToolkit: value
 

test/jumpsystem.jl

@@ -1,4 +1,4 @@
-using ModelingToolkit, DiffEqBase, DiffEqJump, Test, LinearAlgebra
+using ModelingToolkit, DiffEqBase, JumpProcesses, Test, LinearAlgebra
 MT = ModelingToolkit
 
 # basic MT SIR model with tweaks

test/odesystem.jl

@@ -608,18 +608,17 @@ let
     @test sol[y] ≈ 0.9 * sol[x[1]] + sol[x[2]]
     @test isapprox(sol[x[1]][end], 1, atol = 1e-3)
 
-    prob = DAEProblem(sys, [D(y) => 0, D(x[1]) => 0, D(x[2]) => 0], [x[1] => 0.5],
+    prob = DAEProblem(sys, [D(y) => 0, D(x[1]) => 0, D(x[2]) => 0], Pair[x[1] => 0.5],
                       (0, 50))
-    u0_dict = Dict(x[1] => 0.5, x[2] => 0.0)
-    @test prob.u0 ≈ [u0_dict[x] for x in states(sys)]
+    @test prob.u0 ≈ [0.5, 0]
     @test prob.du0 ≈ [0, 0]
     @test prob.p ≈ [1]
     sol = solve(prob, IDA())
     @test isapprox(sol[x[1]][end], 1, atol = 1e-3)
 
     prob = DAEProblem(sys, [D(y) => 0, D(x[1]) => 0, D(x[2]) => 0], Pair[x[1] => 0.5],
                       (0, 50), [k => 2])
-    @test prob.u0 ≈ [u0_dict[x] for x in states(sys)]
+    @test prob.u0 ≈ [0.5, 0]
     @test prob.du0 ≈ [0, 0]
     @test prob.p ≈ [2]
     sol = solve(prob, IDA())