From 8f4061129610f7a0f14f6b0ae779d1a64e71c9ec Mon Sep 17 00:00:00 2001 From: Aayush Sabharwal Date: Mon, 1 Dec 2025 21:12:54 +0530 Subject: [PATCH 01/15] chore!: bump version to 1.0 --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 9683471..e9d862a 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "StateSelection" uuid = "64909d44-ed92-46a8-bbd9-f047dfbdc84b" -version = "0.2.1" authors = ["JuliaHub", "Inc. and other contributors"] +version = "1.0.0" [deps] DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae" From 021b5b07f63ce6b845540d1a8dc155a5d6e2a727 Mon Sep 17 00:00:00 2001 From: Aayush Sabharwal Date: Mon, 1 Dec 2025 21:13:14 +0530 Subject: [PATCH 02/15] refactor: use SciML/BipartiteGraphs.jl --- Project.toml | 14 +- ext/StateSelectionDeepDiffsExt.jl | 6 +- src/StateSelection.jl | 5 +- src/debug.jl | 3 +- src/graph/bipartite.jl | 836 ------------------------------ src/graph/diff.jl | 4 +- src/pantelides.jl | 2 +- src/partial_state_selection.jl | 2 +- src/singularity_removal.jl | 2 +- src/utils.jl | 2 +- 10 files changed, 18 insertions(+), 858 deletions(-) delete mode 100644 src/graph/bipartite.jl diff --git a/Project.toml b/Project.toml index e9d862a..ffe51ad 100644 --- a/Project.toml +++ b/Project.toml @@ -4,6 +4,7 @@ authors = ["JuliaHub", "Inc. and other contributors"] version = "1.0.0" [deps] +BipartiteGraphs = "caf10ac8-0290-4205-88aa-f15908547e8d" DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae" FindFirstFunctions = "64ca27bc-2ba2-4a57-88aa-44e436879224" Graphs = "86223c79-3864-5bf0-83f7-82e725a168b6" @@ -12,7 +13,14 @@ Setfield = "efcf1570-3423-57d1-acb7-fd33fddbac46" SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" UnPack = "3a884ed6-31ef-47d7-9d2a-63182c4928ed" +[weakdeps] +DeepDiffs = "ab62b9b5-e342-54a8-a765-a90f495de1a6" + +[extensions] +StateSelectionDeepDiffsExt = "DeepDiffs" + [compat] +BipartiteGraphs = "0.1.1" DocStringExtensions = "0.9.3" FindFirstFunctions = "1.2.0" Graphs = "1.10.0" @@ -22,12 +30,6 @@ SparseArrays = "1.11.0" UnPack = "1.0.2" julia = "1.9" -[weakdeps] -DeepDiffs = "ab62b9b5-e342-54a8-a765-a90f495de1a6" - -[extensions] -StateSelectionDeepDiffsExt = "DeepDiffs" - [extras] Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" diff --git a/ext/StateSelectionDeepDiffsExt.jl b/ext/StateSelectionDeepDiffsExt.jl index b6b1933..edea1e0 100644 --- a/ext/StateSelectionDeepDiffsExt.jl +++ b/ext/StateSelectionDeepDiffsExt.jl @@ -1,9 +1,7 @@ module StateSelectionDeepDiffsExt using DeepDiffs -using StateSelection.BipartiteGraphs: Label, - BipartiteAdjacencyList, unassigned, - HighlightInt +using BipartiteGraphs: Label, BipartiteAdjacencyList, unassigned, HighlightInt using StateSelection: SystemStructure, MatchedSystemStructure, SystemStructurePrintMatrix @@ -187,4 +185,4 @@ function DeepDiffs.deepdiff(old::Union{MatchedSystemStructure, SystemStructure}, Base.print_matrix(stdout, SystemStructureDiffPrintMatrix(new_sspm, old_sspm)) end -end \ No newline at end of file +end diff --git a/src/StateSelection.jl b/src/StateSelection.jl index 26fee2a..a6cd9ce 100644 --- a/src/StateSelection.jl +++ b/src/StateSelection.jl @@ -6,11 +6,8 @@ using UnPack: @unpack using Graphs # Graph Types -function invview end -function complete end -include("graph/bipartite.jl") +using BipartiteGraphs include("graph/diff.jl") -using .BipartiteGraphs # Math library include("math/bareiss.jl") diff --git a/src/debug.jl b/src/debug.jl index 8e6100c..d9b5a78 100644 --- a/src/debug.jl +++ b/src/debug.jl @@ -1,5 +1,4 @@ - -using .BipartiteGraphs: Label, BipartiteAdjacencyList +using BipartiteGraphs: Label, BipartiteAdjacencyList struct SSAUses{T} eqs::Vector{T} # equation uses diff --git a/src/graph/bipartite.jl b/src/graph/bipartite.jl deleted file mode 100644 index 96dcbee..0000000 --- a/src/graph/bipartite.jl +++ /dev/null @@ -1,836 +0,0 @@ -module BipartiteGraphs - - -export BipartiteEdge, BipartiteGraph, DiCMOBiGraph, Unassigned, unassigned, - Matching, InducedCondensationGraph, maximal_matching, - construct_augmenting_path!, MatchedCondensationGraph - -export 𝑠vertices, 𝑑vertices, has_𝑠vertex, has_𝑑vertex, 𝑠neighbors, 𝑑neighbors, - 𝑠edges, 𝑑edges, nsrcs, ndsts, SRC, DST, set_neighbors!, invview, - delete_srcs!, delete_dsts! - -import ..invview, ..complete - -using DocStringExtensions -using UnPack -using SparseArrays -using Graphs -using Setfield - -### Matching -struct Unassigned - global unassigned - const unassigned = Unassigned.instance -end -# Behaves as a scalar -Base.length(u::Unassigned) = 1 -Base.size(u::Unassigned) = () -Base.iterate(u::Unassigned) = (unassigned, nothing) -Base.iterate(u::Unassigned, state) = nothing - -Base.show(io::IO, ::Unassigned) = printstyled(io, "u"; color = :light_black) - -struct Matching{U, V <: AbstractVector} <: AbstractVector{Union{U, Int}} #=> :Unassigned =# - match::V - inv_match::Union{Nothing, V} -end -# These constructors work around https://github.com/JuliaLang/julia/issues/41948 -Matching{V}(m::Matching) where {V} = convert(Matching{V}, m) -function Base.convert(T::Type{<:Matching{V}}, m::Matching) where {V} - eltype(m) === Union{V, Int} && return M - VUT = typeof(similar(m.match, Union{V, Int})) - Matching{V}(convert(VUT, m.match), - m.inv_match === nothing ? nothing : convert(VUT, m.inv_match)) -end - -Matching(m::Matching) = m -Matching{U}(v::V) where {U, V <: AbstractVector} = Matching{U, V}(v, nothing) -function Matching{U}(v::V, iv::Union{V, Nothing}) where {U, V <: AbstractVector} - Matching{U, V}(v, iv) -end -function Matching(v::V) where {U, V <: AbstractVector{Union{U, Int}}} - Matching{@isdefined(U) ? U : Unassigned, V}(v, nothing) -end -function Matching(m::Int) - Matching{Unassigned}(Union{Int, Unassigned}[unassigned for _ in 1:m], nothing) -end -function Matching{U}(m::Int) where {U} - Matching{Union{Unassigned, U}}(Union{Int, Unassigned, U}[unassigned for _ in 1:m], - nothing) -end - -Base.size(m::Matching) = Base.size(m.match) -Base.getindex(m::Matching, i::Integer) = m.match[i] -Base.iterate(m::Matching, state...) = iterate(m.match, state...) -function Base.copy(m::Matching{U}) where {U} - Matching{U}(copy(m.match), m.inv_match === nothing ? nothing : copy(m.inv_match)) -end -function Base.setindex!(m::Matching{U}, v::Union{Integer, U}, i::Integer) where {U} - if m.inv_match !== nothing - oldv = m.match[i] - # TODO: maybe default Matching to always have an `inv_match`? - - # To maintain the invariant that `m.inv_match[m.match[i]] == i`, we need - # to unassign the matching at `m.inv_match[v]` if it exists. - if v isa Int && (iv = m.inv_match[v]) isa Int - m.match[iv] = unassigned - end - if isa(oldv, Int) - @assert m.inv_match[oldv] == i - m.inv_match[oldv] = unassigned - end - isa(v, Int) && (m.inv_match[v] = i) - end - return m.match[i] = v -end - -function Base.push!(m::Matching, v) - push!(m.match, v) - if v isa Integer && m.inv_match !== nothing - m.inv_match[v] = length(m.match) - end -end - -function complete(m::Matching{U}, - N = maximum((x for x in m.match if isa(x, Int)); init = 0)) where {U} - m.inv_match !== nothing && return m - inv_match = Union{U, Int}[unassigned for _ in 1:N] - for (i, eq) in enumerate(m.match) - isa(eq, Int) || continue - inv_match[eq] = i - end - return Matching{U}(collect(m.match), inv_match) -end - -@noinline function require_complete(m::Matching) - m.inv_match === nothing && - throw(ArgumentError("Backwards matching not defined. `complete` the matching first.")) -end - -function invview(m::Matching{U, V}) where {U, V} - require_complete(m) - return Matching{U, V}(m.inv_match, m.match) -end - -### -### Edges & Vertex -### -@enum VertType SRC DST - -struct BipartiteEdge{I <: Integer} <: Graphs.AbstractEdge{I} - src::I - dst::I - function BipartiteEdge(src::I, dst::V) where {I, V} - T = promote_type(I, V) - new{T}(T(src), T(dst)) - end -end - -Graphs.src(edge::BipartiteEdge) = edge.src -Graphs.dst(edge::BipartiteEdge) = edge.dst - -function Base.show(io::IO, edge::BipartiteEdge) - @unpack src, dst = edge - print(io, "[src: ", src, "] => [dst: ", dst, "]") -end - -Base.:(==)(a::BipartiteEdge, b::BipartiteEdge) = src(a) == src(b) && dst(a) == dst(b) - -### -### Graph -### -""" -$(TYPEDEF) - -A bipartite graph representation between two, possibly distinct, sets of vertices -(source and dependencies). Maps source vertices, labelled `1:N₁`, to vertices -on which they depend (labelled `1:N₂`). - -# Fields -$(FIELDS) - -# Example -```julia -using BipartiteGraphs - -ne = 4 -srcverts = 1:4 -depverts = 1:2 - -# six source vertices -fadjlist = [[1],[1],[2],[2],[1],[1,2]] - -# two vertices they depend on -badjlist = [[1,2,5,6],[3,4,6]] - -bg = BipartiteGraph(7, fadjlist, badjlist) -``` -""" -mutable struct BipartiteGraph{I <: Integer, M} <: Graphs.AbstractGraph{I} - ne::Int - fadjlist::Vector{Vector{I}} # `fadjlist[src] => dsts` - badjlist::Union{Vector{Vector{I}}, I} # `badjlist[dst] => srcs` or `ndsts` - metadata::M -end -function BipartiteGraph(ne::Integer, fadj::AbstractVector, - badj::Union{AbstractVector, Integer} = maximum(maximum, fadj); - metadata = nothing) - BipartiteGraph(ne, fadj, badj, metadata) -end -function BipartiteGraph(fadj::AbstractVector, - badj::Union{AbstractVector, Integer} = maximum(maximum, fadj); - metadata = nothing) - BipartiteGraph(mapreduce(length, +, fadj; init = 0), fadj, badj, metadata) -end - -@noinline function require_complete(g::BipartiteGraph) - g.badjlist isa AbstractVector || - throw(ArgumentError("The graph has no back edges. Use `complete`.")) -end - -function invview(g::BipartiteGraph) - require_complete(g) - BipartiteGraph(g.ne, g.badjlist, g.fadjlist) -end - -function complete(g::BipartiteGraph{I}) where {I} - isa(g.badjlist, AbstractVector) && return g - badjlist = Vector{I}[Vector{I}() for _ in 1:(g.badjlist)] - for (s, l) in enumerate(g.fadjlist) - for d in l - push!(badjlist[d], s) - end - end - BipartiteGraph(g.ne, g.fadjlist, badjlist) -end - -# Matrix whose only purpose is to pretty-print the bipartite graph -struct BipartiteAdjacencyList - u::Union{Vector{Int}, Nothing} - highlight_u::Union{Set{Int}, Nothing} - match # Int or Union{SelectedState, Unassigned} -end -function BipartiteAdjacencyList(u::Union{Vector{Int}, Nothing}) - BipartiteAdjacencyList(u, nothing, unassigned) -end - -function overview_label end -overview_label(::Type) = error("No label defined for this matching type") -overview_label(x) = overview_label(typeof(x)) - -struct HighlightInt - i::Int - highlight::Symbol - match::Bool -end -Base.typeinfo_implicit(::Type{HighlightInt}) = true -function Base.show(io::IO, hi::HighlightInt) - if hi.match - printstyled(io, "(", color = hi.highlight) - printstyled(io, hi.i, color = hi.highlight) - printstyled(io, ")", color = hi.highlight) - else - printstyled(io, hi.i, color = hi.highlight) - end -end - -function Base.show(io::IO, l::BipartiteAdjacencyList) - if !isa(l.match, Union{Int, Unassigned}) - (label, _, color) = overview_label(l.match) - printstyled(io, string(label, " "); color) - else - printstyled(io, " ") - end - if l.u === nothing - printstyled(io, '⋅', color = :light_black) - elseif isempty(l.u) - printstyled(io, '∅', color = :light_black) - elseif l.highlight_u === nothing - print(io, l.u) - else - match = l.match - !isa(match, Int) && (match = unassigned) - function choose_color(i) - solvable = i in l.highlight_u - matched = i == match - if !matched && solvable - :default - elseif !matched && !solvable - :light_black - elseif matched && solvable - :light_yellow - elseif matched && !solvable - :magenta - end - end - if !isempty(setdiff(l.highlight_u, l.u)) - # Only for debugging, shouldn't happen in practice - print(io, - map(union(l.u, l.highlight_u)) do i - HighlightInt(i, !(i in l.u) ? :light_red : choose_color(i), - i == match) - end) - else - print(io, map(l.u) do i - HighlightInt(i, choose_color(i), i == match) - end) - end - end -end - -struct Label - s::String - c::Symbol -end -Label(s::AbstractString) = Label(s, :nothing) -Label(x::Integer) = Label(string(x)) -Base.show(io::IO, l::Label) = printstyled(io, l.s, color = l.c) - -struct BipartiteGraphPrintMatrix <: - AbstractMatrix{Union{Label, Int, BipartiteAdjacencyList}} - bpg::BipartiteGraph -end -Base.size(bgpm::BipartiteGraphPrintMatrix) = (max(nsrcs(bgpm.bpg), ndsts(bgpm.bpg)) + 1, 3) -function Base.getindex(bgpm::BipartiteGraphPrintMatrix, i::Integer, j::Integer) - checkbounds(bgpm, i, j) - if i == 1 - return (Label.(("#", "src", "dst")))[j] - elseif j == 1 - return i - 1 - elseif j == 2 - return BipartiteAdjacencyList(i - 1 <= nsrcs(bgpm.bpg) ? - 𝑠neighbors(bgpm.bpg, i - 1) : nothing) - elseif j == 3 - return BipartiteAdjacencyList(i - 1 <= ndsts(bgpm.bpg) ? - 𝑑neighbors(bgpm.bpg, i - 1) : nothing) - else - @assert false - end -end - -function Base.show(io::IO, b::BipartiteGraph) - print(io, "BipartiteGraph with (", length(b.fadjlist), ", ", - isa(b.badjlist, Int) ? b.badjlist : length(b.badjlist), ") (𝑠,𝑑)-vertices\n") - Base.print_matrix(io, BipartiteGraphPrintMatrix(b)) -end - -""" -```julia -Base.isequal(bg1::BipartiteGraph{T}, bg2::BipartiteGraph{T}) where {T <: Integer} -``` - -Test whether two [`BipartiteGraph`](@ref)s are equal. -""" -function Base.isequal(bg1::BipartiteGraph{T}, bg2::BipartiteGraph{T}) where {T <: Integer} - iseq = (bg1.ne == bg2.ne) - iseq &= (bg1.fadjlist == bg2.fadjlist) - iseq &= (bg1.badjlist == bg2.badjlist) - iseq -end - -""" -$(SIGNATURES) - -Build an empty `BipartiteGraph` with `nsrcs` sources and `ndsts` destinations. -""" -function BipartiteGraph(nsrcs::T, ndsts::T, backedge::Val{B} = Val(true); - metadata = nothing) where {T, B} - fadjlist = map(_ -> T[], 1:nsrcs) - badjlist = B ? map(_ -> T[], 1:ndsts) : ndsts - BipartiteGraph(0, fadjlist, badjlist, metadata) -end - -function Base.copy(bg::BipartiteGraph) - BipartiteGraph(bg.ne, map(copy, bg.fadjlist), map(copy, bg.badjlist), - deepcopy(bg.metadata)) -end -Base.eltype(::Type{<:BipartiteGraph{I}}) where {I} = I -function Base.empty!(g::BipartiteGraph) - foreach(empty!, g.fadjlist) - g.badjlist isa AbstractVector && foreach(empty!, g.badjlist) - g.ne = 0 - if g.metadata !== nothing - foreach(empty!, g.metadata) - end - g -end -Base.length(::BipartiteGraph) = error("length is not well defined! Use `ne` or `nv`.") - -if isdefined(Graphs, :has_contiguous_vertices) - Graphs.has_contiguous_vertices(::Type{<:BipartiteGraph}) = false -end -Graphs.is_directed(::Type{<:BipartiteGraph}) = false -Graphs.vertices(g::BipartiteGraph) = (𝑠vertices(g), 𝑑vertices(g)) -𝑠vertices(g::BipartiteGraph) = axes(g.fadjlist, 1) -function 𝑑vertices(g::BipartiteGraph) - g.badjlist isa AbstractVector ? axes(g.badjlist, 1) : Base.OneTo(g.badjlist) -end -has_𝑠vertex(g::BipartiteGraph, v::Integer) = v in 𝑠vertices(g) -has_𝑑vertex(g::BipartiteGraph, v::Integer) = v in 𝑑vertices(g) -function 𝑠neighbors(g::BipartiteGraph, i::Integer, - with_metadata::Val{M} = Val(false)) where {M} - M ? zip(g.fadjlist[i], g.metadata[i]) : g.fadjlist[i] -end -function 𝑑neighbors(g::BipartiteGraph, j::Integer, - with_metadata::Val{M} = Val(false)) where {M} - require_complete(g) - M ? zip(g.badjlist[j], (g.metadata[i][j] for i in g.badjlist[j])) : g.badjlist[j] -end -Graphs.ne(g::BipartiteGraph) = g.ne -Graphs.nv(g::BipartiteGraph) = sum(length, vertices(g)) -Graphs.edgetype(g::BipartiteGraph{I}) where {I} = BipartiteEdge{I} - -nsrcs(g::BipartiteGraph) = length(𝑠vertices(g)) -ndsts(g::BipartiteGraph) = length(𝑑vertices(g)) - -function Graphs.has_edge(g::BipartiteGraph, edge::BipartiteEdge) - @unpack src, dst = edge - (src in 𝑠vertices(g) && dst in 𝑑vertices(g)) || return false # edge out of bounds - insorted(dst, 𝑠neighbors(g, src)) -end -Base.in(edge::BipartiteEdge, g::BipartiteGraph) = Graphs.has_edge(g, edge) - -### Maximal matching -""" - construct_augmenting_path!(m::Matching, g::BipartiteGraph, vsrc, dstfilter, vcolor=falses(ndsts(g)), ecolor=nothing) -> path_found::Bool - -Try to construct an augmenting path in matching and if such a path is found, -update the matching accordingly. -""" -function construct_augmenting_path!(matching::Matching, g::BipartiteGraph, vsrc, dstfilter, - dcolor = falses(ndsts(g)), scolor = nothing) - scolor === nothing || (scolor[vsrc] = true) - - # if a `vdst` is unassigned and the edge `vsrc <=> vdst` exists - for vdst in 𝑠neighbors(g, vsrc) - if dstfilter(vdst) && matching[vdst] === unassigned - matching[vdst] = vsrc - return true - end - end - - # for every `vsrc` such that edge `vsrc <=> vdst` exists and `vdst` is uncolored - for vdst in 𝑠neighbors(g, vsrc) - (dstfilter(vdst) && !dcolor[vdst]) || continue - dcolor[vdst] = true - if construct_augmenting_path!(matching, g, matching[vdst], dstfilter, dcolor, - scolor) - matching[vdst] = vsrc - return true - end - end - return false -end - -""" - maximal_matching(g::BipartiteGraph; [srcfilter], [dstfilter]) - -For a bipartite graph `g`, construct a maximal matching of destination to source -vertices, subject to the constraint that vertices for which `srcfilter` or `dstfilter`, -return `false` may not be matched. -""" -function maximal_matching(g::BipartiteGraph, ::Type{U} = Unassigned; - srcfilter = vsrc -> true, - dstfilter = vdst -> true) where {U} - matching = Matching{U}(max(nsrcs(g), ndsts(g))) - foreach(Iterators.filter(srcfilter, 𝑠vertices(g))) do vsrc - construct_augmenting_path!(matching, g, vsrc, dstfilter) - end - return matching -end - -### -### Populate -### -struct NoMetadata end -const NO_METADATA = NoMetadata() - -function Graphs.add_edge!(g::BipartiteGraph, i::Integer, j::Integer, md = NO_METADATA) - add_edge!(g, BipartiteEdge(i, j), md) -end -function Graphs.add_edge!(g::BipartiteGraph, edge::BipartiteEdge, md = NO_METADATA) - @unpack fadjlist, badjlist = g - s, d = src(edge), dst(edge) - (has_𝑠vertex(g, s) && has_𝑑vertex(g, d)) || error("edge ($edge) out of range.") - @inbounds list = fadjlist[s] - index = searchsortedfirst(list, d) - @inbounds (index <= length(list) && list[index] == d) && return false # edge already in graph - insert!(list, index, d) - if md !== NO_METADATA - insert!(g.metadata[s], index, md) - end - - g.ne += 1 - if badjlist isa AbstractVector - @inbounds list = badjlist[d] - index = searchsortedfirst(list, s) - insert!(list, index, s) - end - return true # edge successfully added -end - -function Graphs.rem_edge!(g::BipartiteGraph, i::Integer, j::Integer) - Graphs.rem_edge!(g, BipartiteEdge(i, j)) -end -function Graphs.rem_edge!(g::BipartiteGraph, edge::BipartiteEdge) - @unpack fadjlist, badjlist = g - s, d = src(edge), dst(edge) - (has_𝑠vertex(g, s) && has_𝑑vertex(g, d)) || error("edge ($edge) out of range.") - @inbounds list = fadjlist[s] - index = searchsortedfirst(list, d) - @inbounds (index <= length(list) && list[index] == d) || - error("graph does not have edge $edge") - deleteat!(list, index) - g.ne -= 1 - if badjlist isa AbstractVector - @inbounds list = badjlist[d] - index = searchsortedfirst(list, s) - deleteat!(list, index) - end - return true # edge successfully deleted -end - -function Graphs.add_vertex!(g::BipartiteGraph{T}, type::VertType) where {T} - if type === DST - if g.badjlist isa AbstractVector - push!(g.badjlist, T[]) - return length(g.badjlist) - else - g.badjlist += 1 - return g.badjlist - end - elseif type === SRC - push!(g.fadjlist, T[]) - return length(g.fadjlist) - else - error("type ($type) must be either `DST` or `SRC`") - end -end - -function set_neighbors!(g::BipartiteGraph, i::Integer, new_neighbors) - old_neighbors = g.fadjlist[i] - old_nneighbors = length(old_neighbors) - new_nneighbors = length(new_neighbors) - g.ne += new_nneighbors - old_nneighbors - if isa(g.badjlist, AbstractVector) - for n in old_neighbors - @inbounds list = g.badjlist[n] - index = searchsortedfirst(list, i) - if 1 <= index <= length(list) && list[index] == i - deleteat!(list, index) - end - end - for n in new_neighbors - @inbounds list = g.badjlist[n] - index = searchsortedfirst(list, i) - if !(1 <= index <= length(list) && list[index] == i) - insert!(list, index, i) - end - end - end - if iszero(new_nneighbors) # this handles Tuple as well - # Warning: Aliases old_neighbors - empty!(g.fadjlist[i]) - else - g.fadjlist[i] = unique!(sort(new_neighbors)) - end -end - -function delete_srcs!(g::BipartiteGraph, srcs) - for s in srcs - set_neighbors!(g, s, ()) - end - g -end -delete_dsts!(g::BipartiteGraph, srcs) = delete_srcs!(invview(g), srcs) - -### -### Edges iteration -### -Graphs.edges(g::BipartiteGraph) = BipartiteEdgeIter(g, Val(SRC)) -𝑠edges(g::BipartiteGraph) = BipartiteEdgeIter(g, Val(SRC)) -𝑑edges(g::BipartiteGraph) = BipartiteEdgeIter(g, Val(DST)) - -struct BipartiteEdgeIter{T, G} <: Graphs.AbstractEdgeIter - g::G - type::Val{T} -end - -Base.length(it::BipartiteEdgeIter) = ne(it.g) -Base.eltype(it::BipartiteEdgeIter) = edgetype(it.g) - -function Base.iterate(it::BipartiteEdgeIter{SRC, <:BipartiteGraph{T}}, - state = (1, 1, SRC)) where {T} - @unpack g = it - neqs = nsrcs(g) - neqs == 0 && return nothing - eq, jvar = state - - while eq <= neqs - eq′ = eq - vars = 𝑠neighbors(g, eq′) - if jvar > length(vars) - eq += 1 - jvar = 1 - continue - end - edge = BipartiteEdge(eq′, vars[jvar]) - state = (eq, jvar + 1, SRC) - return edge, state - end - return nothing -end - -function Base.iterate(it::BipartiteEdgeIter{DST, <:BipartiteGraph{T}}, - state = (1, 1, DST)) where {T} - @unpack g = it - nvars = ndsts(g) - nvars == 0 && return nothing - ieq, jvar = state - - while jvar <= nvars - eqs = 𝑑neighbors(g, jvar) - if ieq > length(eqs) - ieq = 1 - jvar += 1 - continue - end - edge = BipartiteEdge(eqs[ieq], jvar) - state = (ieq + 1, jvar, DST) - return edge, state - end - return nothing -end - -### -### Utils -### -function Graphs.incidence_matrix(g::BipartiteGraph, val = true) - I = Int[] - J = Int[] - for i in 𝑠vertices(g), n in 𝑠neighbors(g, i) - push!(I, i) - push!(J, n) - end - S = sparse(I, J, val, nsrcs(g), ndsts(g)) -end - -""" - struct DiCMOBiGraph - -This data structure implements a "directed, contracted, matching-oriented" view of an -original (undirected) bipartite graph. It has two modes, depending on the `Transposed` -flag, which switches the direction of the induced matching. - -Essentially the graph adapter performs two largely orthogonal functions -[`Transposed == true` differences are indicated in square brackets]: - - 1. It pairs an undirected bipartite graph with a matching of the destination vertex. - - This matching is used to induce an orientation on the otherwise undirected graph: - Matched edges pass from destination to source [source to destination], all other edges - pass in the opposite direction. - - 2. It exposes the graph view obtained by contracting the destination [source] vertices - along the matched edges. - -The result of this operation is an induced, directed graph on the source [destination] vertices. -The resulting graph has a few desirable properties. In particular, this graph -is acyclic if and only if the induced directed graph on the original bipartite -graph is acyclic. - -# Hypergraph interpretation - -Consider the bipartite graph `B` as the incidence graph of some hypergraph `H`. -Note that a matching `M` on `B` in the above sense is equivalent to determining -an (1,n)-orientation on the hypergraph (i.e. each directed hyperedge has exactly -one head, but any arbitrary number of tails). In this setting, this is simply -the graph formed by expanding each directed hyperedge into `n` ordinary edges -between the same vertices. -""" -mutable struct DiCMOBiGraph{Transposed, I, G <: BipartiteGraph{I}, M <: Matching} <: - Graphs.AbstractGraph{I} - graph::G - ne::Union{Missing, Int} - matching::M - function DiCMOBiGraph{Transposed}(g::G, ne::Union{Missing, Int}, - m::M) where {Transposed, I, G <: BipartiteGraph{I}, M} - new{Transposed, I, G, M}(g, ne, m) - end -end -function DiCMOBiGraph{Transposed}(g::BipartiteGraph) where {Transposed} - DiCMOBiGraph{Transposed}(g, 0, Matching(ndsts(g))) -end -function DiCMOBiGraph{Transposed}(g::BipartiteGraph, m::M) where {Transposed, M} - DiCMOBiGraph{Transposed}(g, missing, m) -end - -function invview(g::DiCMOBiGraph{Transposed}) where {Transposed} - DiCMOBiGraph{!Transposed}(invview(g.graph), g.ne, invview(g.matching)) -end - -Graphs.is_directed(::Type{<:DiCMOBiGraph}) = true -function Graphs.nv(g::DiCMOBiGraph{Transposed}) where {Transposed} - Transposed ? ndsts(g.graph) : nsrcs(g.graph) -end -function Graphs.vertices(g::DiCMOBiGraph{Transposed}) where {Transposed} - Transposed ? 𝑑vertices(g.graph) : 𝑠vertices(g.graph) -end - -struct CMONeighbors{Transposed, V} - g::DiCMOBiGraph{Transposed} - v::V - function CMONeighbors{Transposed}(g::DiCMOBiGraph{Transposed}, - v::V) where {Transposed, V} - new{Transposed, V}(g, v) - end -end - -Graphs.outneighbors(g::DiCMOBiGraph{false}, v) = CMONeighbors{false}(g, v) -Graphs.inneighbors(g::DiCMOBiGraph{false}, v) = inneighbors(invview(g), v) -Base.iterate(c::CMONeighbors{false}) = iterate(c, (c.g.graph.fadjlist[c.v],)) -function Base.iterate(c::CMONeighbors{false}, (l, state...)) - while true - r = iterate(l, state...) - r === nothing && return nothing - # If this is a matched edge, skip it, it's reversed in the induced - # directed graph. Otherwise, if there is no matching for this destination - # edge, also skip it, since it got deleted in the contraction. - vsrc = c.g.matching[r[1]] - if vsrc === c.v || !isa(vsrc, Int) - state = (r[2],) - continue - end - return vsrc, (l, r[2]) - end -end -Base.length(c::CMONeighbors{false}) = count(_ -> true, c) - -liftint(f, x) = (!isa(x, Int)) ? nothing : f(x) -liftnothing(f, x) = x === nothing ? nothing : f(x) - -_vsrc(c::CMONeighbors{true}) = c.g.matching[c.v] -_neighbors(c::CMONeighbors{true}) = liftint(vsrc -> c.g.graph.fadjlist[vsrc], _vsrc(c)) -Base.length(c::CMONeighbors{true}) = something(liftnothing(length, _neighbors(c)), 1) - 1 -Graphs.inneighbors(g::DiCMOBiGraph{true}, v) = CMONeighbors{true}(g, v) -Graphs.outneighbors(g::DiCMOBiGraph{true}, v) = outneighbors(invview(g), v) -Base.iterate(c::CMONeighbors{true}) = liftnothing(ns -> iterate(c, (ns,)), _neighbors(c)) -function Base.iterate(c::CMONeighbors{true}, (l, state...)) - while true - r = iterate(l, state...) - r === nothing && return nothing - if r[1] === c.v - state = (r[2],) - continue - end - return r[1], (l, r[2]) - end -end - -function _edges(g::DiCMOBiGraph{Transposed}) where {Transposed} - Transposed ? - ((w => v for w in inneighbors(g, v)) for v in vertices(g)) : - ((v => w for w in outneighbors(g, v)) for v in vertices(g)) -end - -Graphs.edges(g::DiCMOBiGraph) = (Graphs.SimpleEdge(p) for p in Iterators.flatten(_edges(g))) -function Graphs.ne(g::DiCMOBiGraph) - if g.ne === missing - g.ne = mapreduce(x -> length(x.iter), +, _edges(g)) - end - return g.ne -end - -Graphs.has_edge(g::DiCMOBiGraph{true}, a, b) = a in inneighbors(g, b) -Graphs.has_edge(g::DiCMOBiGraph{false}, a, b) = b in outneighbors(g, a) -# This definition is required for `induced_subgraph` to work -(::Type{<:DiCMOBiGraph})(n::Integer) = SimpleDiGraph(n) - -# Condensation Graphs -abstract type AbstractCondensationGraph <: AbstractGraph{Int} end -function (T::Type{<:AbstractCondensationGraph})(g, sccs::Vector{Union{Int, Vector{Int}}}) - scc_assignment = Vector{Int}(undef, isa(g, BipartiteGraph) ? ndsts(g) : nv(g)) - for (i, c) in enumerate(sccs) - for v in c - scc_assignment[v] = i - end - end - T(g, sccs, scc_assignment) -end -function (T::Type{<:AbstractCondensationGraph})(g, sccs::Vector{Vector{Int}}) - T(g, Vector{Union{Int, Vector{Int}}}(sccs)) -end - -Graphs.is_directed(::Type{<:AbstractCondensationGraph}) = true -Graphs.nv(icg::AbstractCondensationGraph) = length(icg.sccs) -Graphs.vertices(icg::AbstractCondensationGraph) = Base.OneTo(nv(icg)) - -""" - struct MatchedCondensationGraph - -For some bipartite-graph and an orientation induced on its destination contraction, -records the condensation DAG of the digraph formed by the orientation. I.e. this -is a DAG of connected components formed by the destination vertices of some -underlying bipartite graph. -N.B.: This graph does not store explicit neighbor relations of the sccs. -Therefor, the edge multiplicity is derived from the underlying bipartite graph, -i.e. this graph is not strict. -""" -struct MatchedCondensationGraph{G <: DiCMOBiGraph} <: AbstractCondensationGraph - graph::G - # Records the members of a strongly connected component. For efficiency, - # trivial sccs (with one vertex member) are stored inline. Note: the sccs - # here need not be stored in topological order. - sccs::Vector{Union{Int, Vector{Int}}} - # Maps the vertices back to the scc of which they are a part - scc_assignment::Vector{Int} -end - -function Graphs.outneighbors(mcg::MatchedCondensationGraph, cc::Integer) - Iterators.flatten((mcg.scc_assignment[v′] - for v′ in outneighbors(mcg.graph, v) if mcg.scc_assignment[v′] != cc) - for v in mcg.sccs[cc]) -end - -function Graphs.inneighbors(mcg::MatchedCondensationGraph, cc::Integer) - Iterators.flatten((mcg.scc_assignment[v′] - for v′ in inneighbors(mcg.graph, v) if mcg.scc_assignment[v′] != cc) - for v in mcg.sccs[cc]) -end - -""" - struct InducedCondensationGraph - -For some bipartite-graph and a topologicall sorted list of connected components, -represents the condensation DAG of the digraph formed by the orientation. I.e. this -is a DAG of connected components formed by the destination vertices of some -underlying bipartite graph. -N.B.: This graph does not store explicit neighbor relations of the sccs. -Therefor, the edge multiplicity is derived from the underlying bipartite graph, -i.e. this graph is not strict. -""" -struct InducedCondensationGraph{G <: BipartiteGraph} <: AbstractCondensationGraph - graph::G - # Records the members of a strongly connected component. For efficiency, - # trivial sccs (with one vertex member) are stored inline. Note: the sccs - # here are stored in topological order. - sccs::Vector{Union{Int, Vector{Int}}} - # Maps the vertices back to the scc of which they are a part - scc_assignment::Vector{Int} -end - -function _neighbors(icg::InducedCondensationGraph, cc::Integer) - Iterators.flatten(Iterators.flatten(icg.graph.fadjlist[vsrc] - for vsrc in icg.graph.badjlist[v]) - for v in icg.sccs[cc]) -end - -function Graphs.outneighbors(icg::InducedCondensationGraph, v::Integer) - (icg.scc_assignment[n] for n in _neighbors(icg, v) if icg.scc_assignment[n] > v) -end - -function Graphs.inneighbors(icg::InducedCondensationGraph, v::Integer) - (icg.scc_assignment[n] for n in _neighbors(icg, v) if icg.scc_assignment[n] < v) -end - -end # module diff --git a/src/graph/diff.jl b/src/graph/diff.jl index eb7e8a1..bc5a860 100644 --- a/src/graph/diff.jl +++ b/src/graph/diff.jl @@ -72,7 +72,7 @@ function Base.setindex!(dg::DiffGraph, val::Union{Integer, Nothing}, var::Intege end Base.iterate(dg::DiffGraph, state...) = iterate(dg.primal_to_diff, state...) -function complete(dg::DiffGraph) +function BipartiteGraphs.complete(dg::DiffGraph) dg.diff_to_primal !== nothing && return dg diff_to_primal = Union{Int, Nothing}[nothing for _ in 1:length(dg.primal_to_diff)] for (var, diff) in edges(dg) @@ -81,7 +81,7 @@ function complete(dg::DiffGraph) return DiffGraph(dg.primal_to_diff, diff_to_primal) end -function invview(dg::DiffGraph) +function BipartiteGraphs.invview(dg::DiffGraph) require_complete(dg) return DiffGraph(dg.diff_to_primal, dg.primal_to_diff) end diff --git a/src/pantelides.jl b/src/pantelides.jl index e4b64e7..53e2129 100644 --- a/src/pantelides.jl +++ b/src/pantelides.jl @@ -1,4 +1,4 @@ -using .BipartiteGraphs: 𝑑neighbors, 𝑠neighbors, nsrcs, ndsts, +using BipartiteGraphs: 𝑑neighbors, 𝑠neighbors, nsrcs, ndsts, construct_augmenting_path!, unassigned, DiCMOBiGraph """ diff --git a/src/partial_state_selection.jl b/src/partial_state_selection.jl index 2f0d805..5fc420d 100644 --- a/src/partial_state_selection.jl +++ b/src/partial_state_selection.jl @@ -1,4 +1,4 @@ -using .BipartiteGraphs: Unassigned, maximal_matching +using BipartiteGraphs: Unassigned, maximal_matching function partial_state_selection_graph!(state::TransformationState) var_eq_matching = complete(pantelides!(state)) diff --git a/src/singularity_removal.jl b/src/singularity_removal.jl index 989890f..0b61a09 100644 --- a/src/singularity_removal.jl +++ b/src/singularity_removal.jl @@ -1,5 +1,5 @@ using Graphs.Experimental.Traversals -using .BipartiteGraphs: set_neighbors! +using BipartiteGraphs: set_neighbors! function extreme_var(var_to_diff, v, level = nothing, ::Val{descend} = Val(true); callback = _ -> nothing) where {descend} diff --git a/src/utils.jl b/src/utils.jl index a58abb0..e41d295 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -1,7 +1,7 @@ ### ### Bipartite graph utilities ### -using .BipartiteGraphs: 𝑠vertices, 𝑠neighbors +using BipartiteGraphs: 𝑠vertices, 𝑠neighbors n_concrete_eqs(state::TransformationState) = n_concrete_eqs(state.structure) function n_concrete_eqs(graph::BipartiteGraph) From 74df007ea8b4e971233cd56a6d2aa9470dc55796 Mon Sep 17 00:00:00 2001 From: Aayush Sabharwal Date: Mon, 1 Dec 2025 21:13:20 +0530 Subject: [PATCH 03/15] refactor: remove dead code --- src/StructuralTransformations.jl | 73 -------------------------------- 1 file changed, 73 deletions(-) delete mode 100644 src/StructuralTransformations.jl diff --git a/src/StructuralTransformations.jl b/src/StructuralTransformations.jl deleted file mode 100644 index b9aaca3..0000000 --- a/src/StructuralTransformations.jl +++ /dev/null @@ -1,73 +0,0 @@ -module StructuralTransformations - -using Setfield: @set!, @set -using UnPack: @unpack - -using Symbolics: unwrap, linear_expansion, fast_substitute -using SymbolicUtils -using SymbolicUtils.Code -using SymbolicUtils.Rewriters -using SymbolicUtils: similarterm, istree - -using ModelingToolkit -using ModelingToolkit: ODESystem, AbstractSystem, var_from_nested_derivative, Differential, - unknowns, equations, vars, Symbolic, diff2term, value, - operation, arguments, Sym, Term, simplify, solve_for, - isdiffeq, isdifferential, isirreducible, - empty_substitutions, get_substitutions, - get_tearing_state, get_iv, independent_variables, - has_tearing_state, defaults, InvalidSystemException, - ExtraEquationsSystemException, - ExtraVariablesSystemException, - get_postprocess_fbody, vars!, - IncrementalCycleTracker, add_edge_checked!, topological_sort, - invalidate_cache!, Substitutions, get_or_construct_tearing_state, - filter_kwargs, lower_varname, setio, SparseMatrixCLIL, - get_fullvars, has_equations, observed, - Schedule - -using ModelingToolkit.BipartiteGraphs -import .BipartiteGraphs: invview, complete -import ModelingToolkit: var_derivative!, var_derivative_graph! -using Graphs -using ModelingToolkit: algeqs, EquationsView, - SystemStructure, TransformationState, TearingState, - structural_simplify!, - isdiffvar, isdervar, isalgvar, isdiffeq, algeqs, is_only_discrete, - dervars_range, diffvars_range, algvars_range, - DiffGraph, complete!, - get_fullvars, system_subset - -using ModelingToolkit.DiffEqBase -using ModelingToolkit.StaticArrays -using RuntimeGeneratedFunctions: @RuntimeGeneratedFunction, - RuntimeGeneratedFunctions, - drop_expr - -RuntimeGeneratedFunctions.init(@__MODULE__) - -using SparseArrays - -using SimpleNonlinearSolve - -export tearing, partial_state_selection, dae_index_lowering, check_consistency -export dummy_derivative -export build_torn_function, build_observed_function, ODAEProblem -export sorted_incidence_matrix, - pantelides!, pantelides_reassemble, tearing_reassemble, find_solvables!, - linear_subsys_adjmat! -export tearing_assignments, tearing_substitution -export torn_system_jacobian_sparsity -export full_equations -export but_ordered_incidence, lowest_order_variable_mask, highest_order_variable_mask -export computed_highest_diff_variables - -include("utils.jl") -include("pantelides.jl") -include("bipartite_tearing/modia_tearing.jl") -include("tearing.jl") -include("symbolics_tearing.jl") -include("partial_state_selection.jl") -include("codegen.jl") - -end # module From e33acb22b4558554347cb0bcab2fa4ed20115ec1 Mon Sep 17 00:00:00 2001 From: Aayush Sabharwal Date: Mon, 1 Dec 2025 21:14:55 +0530 Subject: [PATCH 04/15] chore: use AGPLv3 license --- LICENSE | 682 ++++++++++++++++++++++++++++++++++++++++++++++++++++++-- 1 file changed, 661 insertions(+), 21 deletions(-) diff --git a/LICENSE b/LICENSE index 4721f7b..be3f7b2 100644 --- a/LICENSE +++ b/LICENSE @@ -1,21 +1,661 @@ -MIT License - -Copyright (c) JuliaHub, Inc. and other contributors - -Permission is hereby granted, free of charge, to any person obtaining a copy -of this software and associated documentation files (the "Software"), to deal -in the Software without restriction, including without limitation the rights -to use, copy, modify, merge, publish, distribute, sublicense, and/or sell -copies of the Software, and to permit persons to whom the Software is -furnished to do so, subject to the following conditions: - -The above copyright notice and this permission notice shall be included in all -copies or substantial portions of the Software. - -THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR -IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, -FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE -AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER -LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, -OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE -SOFTWARE. + GNU AFFERO GENERAL PUBLIC LICENSE + Version 3, 19 November 2007 + + Copyright (C) 2007 Free Software Foundation, Inc. + Everyone is permitted to copy and distribute verbatim copies + of this license document, but changing it is not allowed. + + Preamble + + The GNU Affero General Public License is a free, copyleft license for +software and other kinds of works, specifically designed to ensure +cooperation with the community in the case of network server software. + + The licenses for most software and other practical works are designed +to take away your freedom to share and change the works. By contrast, +our General Public Licenses are intended to guarantee your freedom to +share and change all versions of a program--to make sure it remains free +software for all its users. + + When we speak of free software, we are referring to freedom, not +price. Our General Public Licenses are designed to make sure that you +have the freedom to distribute copies of free software (and charge for +them if you wish), that you receive source code or can get it if you +want it, that you can change the software or use pieces of it in new +free programs, and that you know you can do these things. + + Developers that use our General Public Licenses protect your rights +with two steps: (1) assert copyright on the software, and (2) offer +you this License which gives you legal permission to copy, distribute +and/or modify the software. + + A secondary benefit of defending all users' freedom is that +improvements made in alternate versions of the program, if they +receive widespread use, become available for other developers to +incorporate. Many developers of free software are heartened and +encouraged by the resulting cooperation. However, in the case of +software used on network servers, this result may fail to come about. +The GNU General Public License permits making a modified version and +letting the public access it on a server without ever releasing its +source code to the public. + + The GNU Affero General Public License is designed specifically to +ensure that, in such cases, the modified source code becomes available +to the community. It requires the operator of a network server to +provide the source code of the modified version running there to the +users of that server. Therefore, public use of a modified version, on +a publicly accessible server, gives the public access to the source +code of the modified version. + + An older license, called the Affero General Public License and +published by Affero, was designed to accomplish similar goals. This is +a different license, not a version of the Affero GPL, but Affero has +released a new version of the Affero GPL which permits relicensing under +this license. + + The precise terms and conditions for copying, distribution and +modification follow. + + TERMS AND CONDITIONS + + 0. Definitions. + + "This License" refers to version 3 of the GNU Affero General Public License. + + "Copyright" also means copyright-like laws that apply to other kinds of +works, such as semiconductor masks. + + "The Program" refers to any copyrightable work licensed under this +License. Each licensee is addressed as "you". "Licensees" and +"recipients" may be individuals or organizations. + + To "modify" a work means to copy from or adapt all or part of the work +in a fashion requiring copyright permission, other than the making of an +exact copy. The resulting work is called a "modified version" of the +earlier work or a work "based on" the earlier work. + + A "covered work" means either the unmodified Program or a work based +on the Program. + + To "propagate" a work means to do anything with it that, without +permission, would make you directly or secondarily liable for +infringement under applicable copyright law, except executing it on a +computer or modifying a private copy. Propagation includes copying, +distribution (with or without modification), making available to the +public, and in some countries other activities as well. + + To "convey" a work means any kind of propagation that enables other +parties to make or receive copies. Mere interaction with a user through +a computer network, with no transfer of a copy, is not conveying. + + An interactive user interface displays "Appropriate Legal Notices" +to the extent that it includes a convenient and prominently visible +feature that (1) displays an appropriate copyright notice, and (2) +tells the user that there is no warranty for the work (except to the +extent that warranties are provided), that licensees may convey the +work under this License, and how to view a copy of this License. If +the interface presents a list of user commands or options, such as a +menu, a prominent item in the list meets this criterion. + + 1. Source Code. + + The "source code" for a work means the preferred form of the work +for making modifications to it. "Object code" means any non-source +form of a work. + + A "Standard Interface" means an interface that either is an official +standard defined by a recognized standards body, or, in the case of +interfaces specified for a particular programming language, one that +is widely used among developers working in that language. + + The "System Libraries" of an executable work include anything, other +than the work as a whole, that (a) is included in the normal form of +packaging a Major Component, but which is not part of that Major +Component, and (b) serves only to enable use of the work with that +Major Component, or to implement a Standard Interface for which an +implementation is available to the public in source code form. A +"Major Component", in this context, means a major essential component +(kernel, window system, and so on) of the specific operating system +(if any) on which the executable work runs, or a compiler used to +produce the work, or an object code interpreter used to run it. + + The "Corresponding Source" for a work in object code form means all +the source code needed to generate, install, and (for an executable +work) run the object code and to modify the work, including scripts to +control those activities. However, it does not include the work's +System Libraries, or general-purpose tools or generally available free +programs which are used unmodified in performing those activities but +which are not part of the work. For example, Corresponding Source +includes interface definition files associated with source files for +the work, and the source code for shared libraries and dynamically +linked subprograms that the work is specifically designed to require, +such as by intimate data communication or control flow between those +subprograms and other parts of the work. + + The Corresponding Source need not include anything that users +can regenerate automatically from other parts of the Corresponding +Source. + + The Corresponding Source for a work in source code form is that +same work. + + 2. Basic Permissions. + + All rights granted under this License are granted for the term of +copyright on the Program, and are irrevocable provided the stated +conditions are met. This License explicitly affirms your unlimited +permission to run the unmodified Program. The output from running a +covered work is covered by this License only if the output, given its +content, constitutes a covered work. This License acknowledges your +rights of fair use or other equivalent, as provided by copyright law. + + You may make, run and propagate covered works that you do not +convey, without conditions so long as your license otherwise remains +in force. You may convey covered works to others for the sole purpose +of having them make modifications exclusively for you, or provide you +with facilities for running those works, provided that you comply with +the terms of this License in conveying all material for which you do +not control copyright. Those thus making or running the covered works +for you must do so exclusively on your behalf, under your direction +and control, on terms that prohibit them from making any copies of +your copyrighted material outside their relationship with you. + + Conveying under any other circumstances is permitted solely under +the conditions stated below. Sublicensing is not allowed; section 10 +makes it unnecessary. + + 3. Protecting Users' Legal Rights From Anti-Circumvention Law. + + No covered work shall be deemed part of an effective technological +measure under any applicable law fulfilling obligations under article +11 of the WIPO copyright treaty adopted on 20 December 1996, or +similar laws prohibiting or restricting circumvention of such +measures. + + When you convey a covered work, you waive any legal power to forbid +circumvention of technological measures to the extent such circumvention +is effected by exercising rights under this License with respect to +the covered work, and you disclaim any intention to limit operation or +modification of the work as a means of enforcing, against the work's +users, your or third parties' legal rights to forbid circumvention of +technological measures. + + 4. Conveying Verbatim Copies. + + You may convey verbatim copies of the Program's source code as you +receive it, in any medium, provided that you conspicuously and +appropriately publish on each copy an appropriate copyright notice; +keep intact all notices stating that this License and any +non-permissive terms added in accord with section 7 apply to the code; +keep intact all notices of the absence of any warranty; and give all +recipients a copy of this License along with the Program. + + You may charge any price or no price for each copy that you convey, +and you may offer support or warranty protection for a fee. + + 5. Conveying Modified Source Versions. + + You may convey a work based on the Program, or the modifications to +produce it from the Program, in the form of source code under the +terms of section 4, provided that you also meet all of these conditions: + + a) The work must carry prominent notices stating that you modified + it, and giving a relevant date. + + b) The work must carry prominent notices stating that it is + released under this License and any conditions added under section + 7. This requirement modifies the requirement in section 4 to + "keep intact all notices". + + c) You must license the entire work, as a whole, under this + License to anyone who comes into possession of a copy. This + License will therefore apply, along with any applicable section 7 + additional terms, to the whole of the work, and all its parts, + regardless of how they are packaged. This License gives no + permission to license the work in any other way, but it does not + invalidate such permission if you have separately received it. + + d) If the work has interactive user interfaces, each must display + Appropriate Legal Notices; however, if the Program has interactive + interfaces that do not display Appropriate Legal Notices, your + work need not make them do so. + + A compilation of a covered work with other separate and independent +works, which are not by their nature extensions of the covered work, +and which are not combined with it such as to form a larger program, +in or on a volume of a storage or distribution medium, is called an +"aggregate" if the compilation and its resulting copyright are not +used to limit the access or legal rights of the compilation's users +beyond what the individual works permit. Inclusion of a covered work +in an aggregate does not cause this License to apply to the other +parts of the aggregate. + + 6. Conveying Non-Source Forms. + + You may convey a covered work in object code form under the terms +of sections 4 and 5, provided that you also convey the +machine-readable Corresponding Source under the terms of this License, +in one of these ways: + + a) Convey the object code in, or embodied in, a physical product + (including a physical distribution medium), accompanied by the + Corresponding Source fixed on a durable physical medium + customarily used for software interchange. + + b) Convey the object code in, or embodied in, a physical product + (including a physical distribution medium), accompanied by a + written offer, valid for at least three years and valid for as + long as you offer spare parts or customer support for that product + model, to give anyone who possesses the object code either (1) a + copy of the Corresponding Source for all the software in the + product that is covered by this License, on a durable physical + medium customarily used for software interchange, for a price no + more than your reasonable cost of physically performing this + conveying of source, or (2) access to copy the + Corresponding Source from a network server at no charge. + + c) Convey individual copies of the object code with a copy of the + written offer to provide the Corresponding Source. This + alternative is allowed only occasionally and noncommercially, and + only if you received the object code with such an offer, in accord + with subsection 6b. + + d) Convey the object code by offering access from a designated + place (gratis or for a charge), and offer equivalent access to the + Corresponding Source in the same way through the same place at no + further charge. You need not require recipients to copy the + Corresponding Source along with the object code. If the place to + copy the object code is a network server, the Corresponding Source + may be on a different server (operated by you or a third party) + that supports equivalent copying facilities, provided you maintain + clear directions next to the object code saying where to find the + Corresponding Source. Regardless of what server hosts the + Corresponding Source, you remain obligated to ensure that it is + available for as long as needed to satisfy these requirements. + + e) Convey the object code using peer-to-peer transmission, provided + you inform other peers where the object code and Corresponding + Source of the work are being offered to the general public at no + charge under subsection 6d. + + A separable portion of the object code, whose source code is excluded +from the Corresponding Source as a System Library, need not be +included in conveying the object code work. + + A "User Product" is either (1) a "consumer product", which means any +tangible personal property which is normally used for personal, family, +or household purposes, or (2) anything designed or sold for incorporation +into a dwelling. In determining whether a product is a consumer product, +doubtful cases shall be resolved in favor of coverage. For a particular +product received by a particular user, "normally used" refers to a +typical or common use of that class of product, regardless of the status +of the particular user or of the way in which the particular user +actually uses, or expects or is expected to use, the product. A product +is a consumer product regardless of whether the product has substantial +commercial, industrial or non-consumer uses, unless such uses represent +the only significant mode of use of the product. + + "Installation Information" for a User Product means any methods, +procedures, authorization keys, or other information required to install +and execute modified versions of a covered work in that User Product from +a modified version of its Corresponding Source. The information must +suffice to ensure that the continued functioning of the modified object +code is in no case prevented or interfered with solely because +modification has been made. + + If you convey an object code work under this section in, or with, or +specifically for use in, a User Product, and the conveying occurs as +part of a transaction in which the right of possession and use of the +User Product is transferred to the recipient in perpetuity or for a +fixed term (regardless of how the transaction is characterized), the +Corresponding Source conveyed under this section must be accompanied +by the Installation Information. But this requirement does not apply +if neither you nor any third party retains the ability to install +modified object code on the User Product (for example, the work has +been installed in ROM). + + The requirement to provide Installation Information does not include a +requirement to continue to provide support service, warranty, or updates +for a work that has been modified or installed by the recipient, or for +the User Product in which it has been modified or installed. Access to a +network may be denied when the modification itself materially and +adversely affects the operation of the network or violates the rules and +protocols for communication across the network. + + Corresponding Source conveyed, and Installation Information provided, +in accord with this section must be in a format that is publicly +documented (and with an implementation available to the public in +source code form), and must require no special password or key for +unpacking, reading or copying. + + 7. Additional Terms. + + "Additional permissions" are terms that supplement the terms of this +License by making exceptions from one or more of its conditions. +Additional permissions that are applicable to the entire Program shall +be treated as though they were included in this License, to the extent +that they are valid under applicable law. If additional permissions +apply only to part of the Program, that part may be used separately +under those permissions, but the entire Program remains governed by +this License without regard to the additional permissions. + + When you convey a copy of a covered work, you may at your option +remove any additional permissions from that copy, or from any part of +it. (Additional permissions may be written to require their own +removal in certain cases when you modify the work.) You may place +additional permissions on material, added by you to a covered work, +for which you have or can give appropriate copyright permission. + + Notwithstanding any other provision of this License, for material you +add to a covered work, you may (if authorized by the copyright holders of +that material) supplement the terms of this License with terms: + + a) Disclaiming warranty or limiting liability differently from the + terms of sections 15 and 16 of this License; or + + b) Requiring preservation of specified reasonable legal notices or + author attributions in that material or in the Appropriate Legal + Notices displayed by works containing it; or + + c) Prohibiting misrepresentation of the origin of that material, or + requiring that modified versions of such material be marked in + reasonable ways as different from the original version; or + + d) Limiting the use for publicity purposes of names of licensors or + authors of the material; or + + e) Declining to grant rights under trademark law for use of some + trade names, trademarks, or service marks; or + + f) Requiring indemnification of licensors and authors of that + material by anyone who conveys the material (or modified versions of + it) with contractual assumptions of liability to the recipient, for + any liability that these contractual assumptions directly impose on + those licensors and authors. + + All other non-permissive additional terms are considered "further +restrictions" within the meaning of section 10. If the Program as you +received it, or any part of it, contains a notice stating that it is +governed by this License along with a term that is a further +restriction, you may remove that term. If a license document contains +a further restriction but permits relicensing or conveying under this +License, you may add to a covered work material governed by the terms +of that license document, provided that the further restriction does +not survive such relicensing or conveying. + + If you add terms to a covered work in accord with this section, you +must place, in the relevant source files, a statement of the +additional terms that apply to those files, or a notice indicating +where to find the applicable terms. + + Additional terms, permissive or non-permissive, may be stated in the +form of a separately written license, or stated as exceptions; +the above requirements apply either way. + + 8. Termination. + + You may not propagate or modify a covered work except as expressly +provided under this License. Any attempt otherwise to propagate or +modify it is void, and will automatically terminate your rights under +this License (including any patent licenses granted under the third +paragraph of section 11). + + However, if you cease all violation of this License, then your +license from a particular copyright holder is reinstated (a) +provisionally, unless and until the copyright holder explicitly and +finally terminates your license, and (b) permanently, if the copyright +holder fails to notify you of the violation by some reasonable means +prior to 60 days after the cessation. + + Moreover, your license from a particular copyright holder is +reinstated permanently if the copyright holder notifies you of the +violation by some reasonable means, this is the first time you have +received notice of violation of this License (for any work) from that +copyright holder, and you cure the violation prior to 30 days after +your receipt of the notice. + + Termination of your rights under this section does not terminate the +licenses of parties who have received copies or rights from you under +this License. If your rights have been terminated and not permanently +reinstated, you do not qualify to receive new licenses for the same +material under section 10. + + 9. Acceptance Not Required for Having Copies. + + You are not required to accept this License in order to receive or +run a copy of the Program. Ancillary propagation of a covered work +occurring solely as a consequence of using peer-to-peer transmission +to receive a copy likewise does not require acceptance. However, +nothing other than this License grants you permission to propagate or +modify any covered work. These actions infringe copyright if you do +not accept this License. Therefore, by modifying or propagating a +covered work, you indicate your acceptance of this License to do so. + + 10. Automatic Licensing of Downstream Recipients. + + Each time you convey a covered work, the recipient automatically +receives a license from the original licensors, to run, modify and +propagate that work, subject to this License. You are not responsible +for enforcing compliance by third parties with this License. + + An "entity transaction" is a transaction transferring control of an +organization, or substantially all assets of one, or subdividing an +organization, or merging organizations. If propagation of a covered +work results from an entity transaction, each party to that +transaction who receives a copy of the work also receives whatever +licenses to the work the party's predecessor in interest had or could +give under the previous paragraph, plus a right to possession of the +Corresponding Source of the work from the predecessor in interest, if +the predecessor has it or can get it with reasonable efforts. + + You may not impose any further restrictions on the exercise of the +rights granted or affirmed under this License. For example, you may +not impose a license fee, royalty, or other charge for exercise of +rights granted under this License, and you may not initiate litigation +(including a cross-claim or counterclaim in a lawsuit) alleging that +any patent claim is infringed by making, using, selling, offering for +sale, or importing the Program or any portion of it. + + 11. Patents. + + A "contributor" is a copyright holder who authorizes use under this +License of the Program or a work on which the Program is based. The +work thus licensed is called the contributor's "contributor version". + + A contributor's "essential patent claims" are all patent claims +owned or controlled by the contributor, whether already acquired or +hereafter acquired, that would be infringed by some manner, permitted +by this License, of making, using, or selling its contributor version, +but do not include claims that would be infringed only as a +consequence of further modification of the contributor version. For +purposes of this definition, "control" includes the right to grant +patent sublicenses in a manner consistent with the requirements of +this License. + + Each contributor grants you a non-exclusive, worldwide, royalty-free +patent license under the contributor's essential patent claims, to +make, use, sell, offer for sale, import and otherwise run, modify and +propagate the contents of its contributor version. + + In the following three paragraphs, a "patent license" is any express +agreement or commitment, however denominated, not to enforce a patent +(such as an express permission to practice a patent or covenant not to +sue for patent infringement). To "grant" such a patent license to a +party means to make such an agreement or commitment not to enforce a +patent against the party. + + If you convey a covered work, knowingly relying on a patent license, +and the Corresponding Source of the work is not available for anyone +to copy, free of charge and under the terms of this License, through a +publicly available network server or other readily accessible means, +then you must either (1) cause the Corresponding Source to be so +available, or (2) arrange to deprive yourself of the benefit of the +patent license for this particular work, or (3) arrange, in a manner +consistent with the requirements of this License, to extend the patent +license to downstream recipients. "Knowingly relying" means you have +actual knowledge that, but for the patent license, your conveying the +covered work in a country, or your recipient's use of the covered work +in a country, would infringe one or more identifiable patents in that +country that you have reason to believe are valid. + + If, pursuant to or in connection with a single transaction or +arrangement, you convey, or propagate by procuring conveyance of, a +covered work, and grant a patent license to some of the parties +receiving the covered work authorizing them to use, propagate, modify +or convey a specific copy of the covered work, then the patent license +you grant is automatically extended to all recipients of the covered +work and works based on it. + + A patent license is "discriminatory" if it does not include within +the scope of its coverage, prohibits the exercise of, or is +conditioned on the non-exercise of one or more of the rights that are +specifically granted under this License. You may not convey a covered +work if you are a party to an arrangement with a third party that is +in the business of distributing software, under which you make payment +to the third party based on the extent of your activity of conveying +the work, and under which the third party grants, to any of the +parties who would receive the covered work from you, a discriminatory +patent license (a) in connection with copies of the covered work +conveyed by you (or copies made from those copies), or (b) primarily +for and in connection with specific products or compilations that +contain the covered work, unless you entered into that arrangement, +or that patent license was granted, prior to 28 March 2007. + + Nothing in this License shall be construed as excluding or limiting +any implied license or other defenses to infringement that may +otherwise be available to you under applicable patent law. + + 12. No Surrender of Others' Freedom. + + If conditions are imposed on you (whether by court order, agreement or +otherwise) that contradict the conditions of this License, they do not +excuse you from the conditions of this License. If you cannot convey a +covered work so as to satisfy simultaneously your obligations under this +License and any other pertinent obligations, then as a consequence you may +not convey it at all. For example, if you agree to terms that obligate you +to collect a royalty for further conveying from those to whom you convey +the Program, the only way you could satisfy both those terms and this +License would be to refrain entirely from conveying the Program. + + 13. Remote Network Interaction; Use with the GNU General Public License. + + Notwithstanding any other provision of this License, if you modify the +Program, your modified version must prominently offer all users +interacting with it remotely through a computer network (if your version +supports such interaction) an opportunity to receive the Corresponding +Source of your version by providing access to the Corresponding Source +from a network server at no charge, through some standard or customary +means of facilitating copying of software. This Corresponding Source +shall include the Corresponding Source for any work covered by version 3 +of the GNU General Public License that is incorporated pursuant to the +following paragraph. + + Notwithstanding any other provision of this License, you have +permission to link or combine any covered work with a work licensed +under version 3 of the GNU General Public License into a single +combined work, and to convey the resulting work. The terms of this +License will continue to apply to the part which is the covered work, +but the work with which it is combined will remain governed by version +3 of the GNU General Public License. + + 14. Revised Versions of this License. + + The Free Software Foundation may publish revised and/or new versions of +the GNU Affero General Public License from time to time. Such new versions +will be similar in spirit to the present version, but may differ in detail to +address new problems or concerns. + + Each version is given a distinguishing version number. If the +Program specifies that a certain numbered version of the GNU Affero General +Public License "or any later version" applies to it, you have the +option of following the terms and conditions either of that numbered +version or of any later version published by the Free Software +Foundation. If the Program does not specify a version number of the +GNU Affero General Public License, you may choose any version ever published +by the Free Software Foundation. + + If the Program specifies that a proxy can decide which future +versions of the GNU Affero General Public License can be used, that proxy's +public statement of acceptance of a version permanently authorizes you +to choose that version for the Program. + + Later license versions may give you additional or different +permissions. However, no additional obligations are imposed on any +author or copyright holder as a result of your choosing to follow a +later version. + + 15. Disclaimer of Warranty. + + THERE IS NO WARRANTY FOR THE PROGRAM, TO THE EXTENT PERMITTED BY +APPLICABLE LAW. EXCEPT WHEN OTHERWISE STATED IN WRITING THE COPYRIGHT +HOLDERS AND/OR OTHER PARTIES PROVIDE THE PROGRAM "AS IS" WITHOUT WARRANTY +OF ANY KIND, EITHER EXPRESSED OR IMPLIED, INCLUDING, BUT NOT LIMITED TO, +THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR +PURPOSE. THE ENTIRE RISK AS TO THE QUALITY AND PERFORMANCE OF THE PROGRAM +IS WITH YOU. SHOULD THE PROGRAM PROVE DEFECTIVE, YOU ASSUME THE COST OF +ALL NECESSARY SERVICING, REPAIR OR CORRECTION. + + 16. Limitation of Liability. + + IN NO EVENT UNLESS REQUIRED BY APPLICABLE LAW OR AGREED TO IN WRITING +WILL ANY COPYRIGHT HOLDER, OR ANY OTHER PARTY WHO MODIFIES AND/OR CONVEYS +THE PROGRAM AS PERMITTED ABOVE, BE LIABLE TO YOU FOR DAMAGES, INCLUDING ANY +GENERAL, SPECIAL, INCIDENTAL OR CONSEQUENTIAL DAMAGES ARISING OUT OF THE +USE OR INABILITY TO USE THE PROGRAM (INCLUDING BUT NOT LIMITED TO LOSS OF +DATA OR DATA BEING RENDERED INACCURATE OR LOSSES SUSTAINED BY YOU OR THIRD +PARTIES OR A FAILURE OF THE PROGRAM TO OPERATE WITH ANY OTHER PROGRAMS), +EVEN IF SUCH HOLDER OR OTHER PARTY HAS BEEN ADVISED OF THE POSSIBILITY OF +SUCH DAMAGES. + + 17. Interpretation of Sections 15 and 16. + + If the disclaimer of warranty and limitation of liability provided +above cannot be given local legal effect according to their terms, +reviewing courts shall apply local law that most closely approximates +an absolute waiver of all civil liability in connection with the +Program, unless a warranty or assumption of liability accompanies a +copy of the Program in return for a fee. + + END OF TERMS AND CONDITIONS + + How to Apply These Terms to Your New Programs + + If you develop a new program, and you want it to be of the greatest +possible use to the public, the best way to achieve this is to make it +free software which everyone can redistribute and change under these terms. + + To do so, attach the following notices to the program. It is safest +to attach them to the start of each source file to most effectively +state the exclusion of warranty; and each file should have at least +the "copyright" line and a pointer to where the full notice is found. + + + Copyright (C) + + This program is free software: you can redistribute it and/or modify + it under the terms of the GNU Affero General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + This program is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU Affero General Public License for more details. + + You should have received a copy of the GNU Affero General Public License + along with this program. If not, see . + +Also add information on how to contact you by electronic and paper mail. + + If your software can interact with users remotely through a computer +network, you should also make sure that it provides a way for users to +get its source. For example, if your program is a web application, its +interface could display a "Source" link that leads users to an archive +of the code. There are many ways you could offer source, and different +solutions will be better for different programs; see section 13 for the +specific requirements. + + You should also get your employer (if you work as a programmer) or school, +if any, to sign a "copyright disclaimer" for the program, if necessary. +For more information on this, and how to apply and follow the GNU AGPL, see +. From 0f3b49bf71421d6d7e1c126efe4073dd3d4c5768 Mon Sep 17 00:00:00 2001 From: Aayush Sabharwal Date: Tue, 2 Dec 2025 16:30:21 +0530 Subject: [PATCH 05/15] feat: improve and document interface functions --- src/interface.jl | 118 ++++++++++++++++++++++++++++++++++++++++++++++- 1 file changed, 117 insertions(+), 1 deletion(-) diff --git a/src/interface.jl b/src/interface.jl index a572fca..46e8f3a 100644 --- a/src/interface.jl +++ b/src/interface.jl @@ -1,12 +1,87 @@ +""" + $TYPEDEF + +Supertype for a mutable struct representing structural information about a DAE. Must have +the following fields: + +- `var_to_diff::DiffGraph`: A mapping from (indices of) variables to (the indices of) their + derivatives. +- `eq_to_diff::DiffGraph`: A mapping from (indices of) equations to (the indices of) their + derivatives. +- `graph::BipartiteGraph{Int, Nothing}`: The bipartite incidence graph of the system. + Source vertices are equations, destination vertices are variables. +- `solvable_graph::Union{BipartiteGraph{Int, Nothing}, Nothing}`: Similar to `graph`, but + instead of incidence it tracks which equations are linearly solvable for which variables. + +Any additional fields are left up to the implementor. +""" abstract type SystemStructure; end is_only_discrete(::SystemStructure) = false +""" + $TYPEDEF + +Supertype for structs representing the state of a system of DAEs during structural +transformations. Must have the following fields: + +- `structure<:SystemStructure`: A `SystemStructure` subtype. Should ideally be + concretely typed to a specific implementation. +- `fullvars::AbstractVector{T}`: A list of variables in the system, ordered identically + to the destination vertices of `structure.graph`. The `eltype` of this buffer can be + chosen by the implementor. If this field is not present, `get_fullvars` must be + implemented for the type. + +In addition to the structural information in `structure`, this struct typically contains +information relevant to the symbolic structure of the system. This can be used to +reconstruct the system for code-generation after structural transformations. +""" abstract type TransformationState{T} end -abstract type AbstractTearingState{T} <: TransformationState{T} end + +""" + $TYPEDSIGNATURES + +Get the ordered list of variables in the given `TransformationState`. Defaults to +`ts.fullvars`. +""" +@inline get_fullvars(ts::TransformationState) = ts.fullvars + +""" + $TYPEDSIGNATURES + +Populate `state.structure.solvable_graph` with information about the solvability of +equations. The implementation should validate that `state.structure.solvable_graph` +is `nothing` before modifying `state`. The default implementation relies on +[`find_eq_solvables!`](@ref), to which it forwards all keyword arguments. +""" +function find_solvables!(state::TransformationState; kwargs...) + @assert state.structure.solvable_graph === nothing + graph = state.structure.graph + state.structure.solvable_graph = BipartiteGraph(nsrcs(graph), ndsts(graph)) + for ieq in 1:nsrcs(graph) + find_eq_solvables!(state, ieq; kwargs...) + end + return nothing +end + +""" + find_eq_solvables!(state::TransformationState, ieq::Int; kwargs...) + +Identify which variables equation `ieq` can be rearranged to solve for, and populate +`state.structure.solvable_graph` accordingly. Keyword arguments are left to the +implementor, and can influence the criteria for solvability. +""" +function find_eq_solvables! end + struct SelectedState end BipartiteGraphs.overview_label(::Type{SelectedState}) = ('∫', " Selected State", :cyan) +""" + linear_subsys_adjmat!(state::TransformationState; kwargs...) + +Find the adjacency matrix of linear subsystems in `state` and return them as a +`SparseMatrixCLIL`. May mutate `state` to cache linearity information. +""" function linear_subsys_adjmat! end function eq_derivative! end function var_derivative! end @@ -29,6 +104,13 @@ function var_derivative_graph!(s::SystemStructure, v::Int) return var_diff end +""" + $TYPEDSIGNATURES + +Call `BipartiteGraphs.complete` on the documented required fields of a `SystemStructure`. +This method may also be implemented manually to perform additional tasks. Should return +the completed `SystemStructure`. +""" function complete!(s::SystemStructure) s.var_to_diff = complete(s.var_to_diff) s.eq_to_diff = complete(s.eq_to_diff) @@ -38,3 +120,37 @@ function complete!(s::SystemStructure) end s end + +isdervar(s::SystemStructure, i) = invview(s.var_to_diff)[i] !== nothing +function isalgvar(s::SystemStructure, i) + s.var_to_diff[i] === nothing && invview(s.var_to_diff)[i] === nothing +end +function isdiffvar(s::SystemStructure, i) + s.var_to_diff[i] !== nothing && invview(s.var_to_diff)[i] === nothing +end + +function dervars_range(s::SystemStructure) + Iterators.filter(Base.Fix1(isdervar, s), Base.OneTo(ndsts(s.graph))) +end +function diffvars_range(s::SystemStructure) + Iterators.filter(Base.Fix1(isdiffvar, s), Base.OneTo(ndsts(s.graph))) +end +function algvars_range(s::SystemStructure) + Iterators.filter(Base.Fix1(isalgvar, s), Base.OneTo(ndsts(s.graph))) +end + +function algeqs(s::SystemStructure) + BitSet(findall(map(1:nsrcs(s.graph)) do eq + all(v -> !isdervar(s, v), 𝑠neighbors(s.graph, eq)) + end)) +end + +""" + $TYPEDSIGNATURES + +Find an equation-variable maximal bipartite matching for `s`, using the incidence graph +`s.graph`. +""" +function BipartiteGraphs.maximal_matching(s::SystemStructure; kw...) + maximal_matching(s.graph; kw...) +end From 8ef1ed00b278416f54375d179364bb7a4f44e363 Mon Sep 17 00:00:00 2001 From: Aayush Sabharwal Date: Tue, 2 Dec 2025 16:33:27 +0530 Subject: [PATCH 06/15] feat: handle non-fully-determined systems in `tear_graph_modia` --- src/modia_tearing.jl | 14 ++++++++++---- 1 file changed, 10 insertions(+), 4 deletions(-) diff --git a/src/modia_tearing.jl b/src/modia_tearing.jl index da29b07..a3c346f 100644 --- a/src/modia_tearing.jl +++ b/src/modia_tearing.jl @@ -91,6 +91,8 @@ function tear_graph_modia(structure::SystemStructure, isder::F = nothing, ieqs = Int[] filtered_vars = BitSet() + free_eqs = free_equations(graph, var_sccs, var_eq_matching, varfilter) + is_overdetemined = !isempty(free_eqs) for vars in var_sccs for var in vars if varfilter(var) @@ -105,10 +107,14 @@ function tear_graph_modia(structure::SystemStructure, isder::F = nothing, filtered_vars, isder) - # clear cache - vargraph.ne = 0 - for var in vars - vargraph.matching[var] = unassigned + # If the systems is overdetemined, we cannot assume the free equations + # will not form algebraic loops with equations in the sccs. + if !is_overdetemined + # clear cache + vargraph.ne = 0 + for var in vars + vargraph.matching[var] = unassigned + end end empty!(ieqs) empty!(filtered_vars) From 0e9c3911eb2aff7eb85c427a96561646d4df6950 Mon Sep 17 00:00:00 2001 From: Aayush Sabharwal Date: Tue, 2 Dec 2025 16:33:50 +0530 Subject: [PATCH 07/15] docs: add docstrings --- src/pantelides.jl | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/src/pantelides.jl b/src/pantelides.jl index 53e2129..0a06196 100644 --- a/src/pantelides.jl +++ b/src/pantelides.jl @@ -2,7 +2,7 @@ using BipartiteGraphs: 𝑑neighbors, 𝑠neighbors, nsrcs, ndsts, construct_augmenting_path!, unassigned, DiCMOBiGraph """ - computed_highest_diff_variables(structure) + $TYPEDSIGNATURES Computes which variables are the "highest-differentiated" for purposes of pantelides. Ordinarily this is relatively straightforward. However, in our @@ -14,8 +14,11 @@ case, there is one complicating condition: This function takes care of these complications are returns a boolean array for every variable, indicating whether it is considered "highest-differentiated". + +`varfilter` is a filter function which takes index of a variable in `structure` and +determines whether it should be included in the list. """ -function computed_highest_diff_variables(structure, varfilter) +function computed_highest_diff_variables(structure::SystemStructure, varfilter) @unpack graph, var_to_diff = structure nvars = length(var_to_diff) From f4ec6b5aae3dc95e2945d472f87edb6b518dc2d3 Mon Sep 17 00:00:00 2001 From: Aayush Sabharwal Date: Tue, 2 Dec 2025 16:34:22 +0530 Subject: [PATCH 08/15] refactor: improve type-stability of `structural_singularity_removal!` --- src/singularity_removal.jl | 16 ++++++++++++---- 1 file changed, 12 insertions(+), 4 deletions(-) diff --git a/src/singularity_removal.jl b/src/singularity_removal.jl index 0b61a09..393241f 100644 --- a/src/singularity_removal.jl +++ b/src/singularity_removal.jl @@ -25,10 +25,12 @@ function structural_singularity_removal!(state::TransformationState; @unpack graph, var_to_diff, solvable_graph = state.structure mm = structural_singularity_removal!(state, mm) s = state.structure - for g in (s.graph, s.solvable_graph) - g === nothing && continue + for (ei, e) in enumerate(mm.nzrows) + set_neighbors!(s.graph, e, mm.row_cols[ei]) + end + if s.solvable_graph isa BipartiteGraph{Int, Nothing} for (ei, e) in enumerate(mm.nzrows) - set_neighbors!(g, e, mm.row_cols[ei]) + set_neighbors!(s.solvable_graph, e, mm.row_cols[ei]) end end @@ -202,7 +204,13 @@ function aag_bareiss!(structure, mm_orig::SparseMatrixCLIL{T, Ti}) where {T, Ti} bar = do_bareiss!(mm, mm_orig, is_linear_variables, is_highest_diff) end - return mm, solvable_variables, bar + # This phrasing infers the return type as `Union{Tuple{...}}` instead of + # `Tuple{Union{...}, ...}` + if mm isa SparseMatrixCLIL{BigInt, Ti} + return mm, solvable_variables, bar + else + return mm, solvable_variables, bar + end end function do_bareiss!(M, Mold, is_linear_variables, is_highest_diff) From 2c2c56c965ce617c4c1bbd8315630cf03b570572 Mon Sep 17 00:00:00 2001 From: Aayush Sabharwal Date: Tue, 2 Dec 2025 16:34:54 +0530 Subject: [PATCH 09/15] feat: document, modularize and add some utility functions --- src/StateSelection.jl | 1 + src/utils.jl | 93 ++++++++++++++++++++++++++++++++++--------- 2 files changed, 76 insertions(+), 18 deletions(-) diff --git a/src/StateSelection.jl b/src/StateSelection.jl index a6cd9ce..3977b7c 100644 --- a/src/StateSelection.jl +++ b/src/StateSelection.jl @@ -4,6 +4,7 @@ using DocStringExtensions using Setfield: @set!, @set using UnPack: @unpack using Graphs +import SparseArrays # Graph Types using BipartiteGraphs diff --git a/src/utils.jl b/src/utils.jl index e41d295..0e6aa6c 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -4,8 +4,9 @@ using BipartiteGraphs: 𝑠vertices, 𝑠neighbors n_concrete_eqs(state::TransformationState) = n_concrete_eqs(state.structure) +n_concrete_eqs(structure::SystemStructure) = n_concrete_eqs(structure.graph) function n_concrete_eqs(graph::BipartiteGraph) - neqs = count(e -> !isempty(𝑠neighbors(graph, e)), 𝑠vertices(graph)) + count(e -> !isempty(𝑠neighbors(graph, e)), 𝑠vertices(graph)) end struct InvalidSystemException <: Exception @@ -75,8 +76,41 @@ end ### ### Structural check ### -function check_consistency(state::TransformationState, orig_inputs) +""" + $(TYPEDSIGNATURES) + +Check if the `state` represents a singular system, and return the unmatched variables. +""" +function singular_check(state::TransformationState) + (; graph, var_to_diff) = state.structure fullvars = get_fullvars(state) + # This is defined to check if Pantelides algorithm terminates. For more + # details, check the equation (15) of the original paper. + extended_graph = (@set graph.fadjlist = Vector{Int}[graph.fadjlist; + map(collect, edges(var_to_diff))]) + extended_var_eq_matching = maximal_matching(extended_graph) + + nvars = ndsts(graph) + unassigned_var = eltype(get_fullvars(state))[] + for (vj, eq) in enumerate(extended_var_eq_matching) + vj > nvars && break + if eq === unassigned && !isempty(𝑑neighbors(graph, vj)) + push!(unassigned_var, fullvars[vj]) + end + end + return unassigned_var + +end + +""" + $(TYPEDSIGNATURES) + +Check the consistency of `state`, given the inputs `orig_inputs`. If `nothrow == false`, +throws an error if the system is under-/over-determined or singular. In this case, if the +function returns it will return `true`. If `nothrow == true`, it will return `false` +instead of throwing an error. The singular case will print a warning. +""" +function check_consistency(state::TransformationState, orig_inputs; nothrow = false) neqs = n_concrete_eqs(state) @unpack graph, var_to_diff = state.structure highest_vars = computed_highest_diff_variables(complete!(state.structure)) @@ -89,6 +123,7 @@ function check_consistency(state::TransformationState, orig_inputs) is_balanced = n_highest_vars == neqs if neqs > 0 && !is_balanced + nothrow && return false varwhitelist = var_to_diff .== nothing var_eq_matching = maximal_matching(graph, dstfilter = v -> varwhitelist[v]) # not assigned @@ -103,20 +138,10 @@ function check_consistency(state::TransformationState, orig_inputs) error_reporting(state, bad_idxs, n_highest_vars, iseqs, orig_inputs) end - # This is defined to check if Pantelides algorithm terminates. For more - # details, check the equation (15) of the original paper. - extended_graph = (@set graph.fadjlist = Vector{Int}[graph.fadjlist; - map(collect, edges(var_to_diff))]) - extended_var_eq_matching = maximal_matching(extended_graph) - - unassigned_var = [] - for (vj, eq) in enumerate(extended_var_eq_matching) - if eq === unassigned && !isempty(𝑑neighbors(graph, vj)) - push!(unassigned_var, fullvars[vj]) - end - end + unassigned_var = singular_check(state) if !isempty(unassigned_var) || !is_balanced + nothrow && return false io = IOBuffer() Base.print_array(io, unassigned_var) unassigned_var_str = String(take!(io)) @@ -126,7 +151,7 @@ function check_consistency(state::TransformationState, orig_inputs) throw(InvalidSystemException(errmsg)) end - return nothing + return true end ### @@ -134,20 +159,52 @@ end ### """ - find_var_sccs(g::BipartiteGraph, assign=nothing) + $TYPEDSIGNATURES Find strongly connected components of the variables defined by `g`. `assign` gives the undirected bipartite graph a direction. When `assign === nothing`, we assume that the ``i``-th variable is assigned to the ``i``-th equation. + +If `topsort == true`, topologically sort the SCCs. """ -function find_var_sccs(g::BipartiteGraph, assign = nothing) +function find_var_sccs(g::BipartiteGraph, assign = nothing; topsort = false) cmog = DiCMOBiGraph{true}(g, Matching(assign === nothing ? Base.OneTo(nsrcs(g)) : assign)) sccs = Graphs.strongly_connected_components(cmog) + if topsort + cgraph = MatchedCondensationGraph(cmog, sccs) + toporder = topological_sort(cgraph) + sccs = sccs[toporder] + end foreach(sort!, sccs) return sccs end +""" + $TYPEDSIGNATURES + +Find strongly connected components of algebraic variables in a system. +""" +function algebraic_variables_scc(structure::SystemStructure) + graph = structure.graph + # skip over differential equations + algvars = BitSet(findall(v -> isalgvar(structure, v), 1:ndsts(graph))) + algeqs = BitSet(findall(map(1:nsrcs(graph)) do eq + all(v -> !isdervar(structure, v), + 𝑠neighbors(graph, eq)) + end)) + var_eq_matching = complete( + maximal_matching(graph, e -> e in algeqs, v -> v in algvars), ndsts(graph)) + var_sccs = find_var_sccs(complete(graph), var_eq_matching) + + return var_eq_matching, var_sccs +end + +""" + $(TYPEDSIGNATURES) + +Obtain the incidence matrix of the system sorted by the algebraic SCCs. +""" function sorted_incidence_matrix(ts::TransformationState, val = true; only_algeqs = false, only_algvars = false) var_eq_matching, var_scc = algebraic_variables_scc(ts) @@ -190,5 +247,5 @@ function sorted_incidence_matrix(ts::TransformationState, val = true; only_algeq push!(J, j) end end - sparse(I, J, val, nsrcs(graph), ndsts(graph)) + SparseArrays.sparse(I, J, val, nsrcs(graph), ndsts(graph)) end From 78951f41693ec6f5901863240a635dead1f54c27 Mon Sep 17 00:00:00 2001 From: Aayush Sabharwal Date: Tue, 2 Dec 2025 16:36:23 +0530 Subject: [PATCH 10/15] feat: add `TearingAlgorithm` interface --- src/tearing.jl | 68 ++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 68 insertions(+) diff --git a/src/tearing.jl b/src/tearing.jl index 54b1ec6..02c1f15 100644 --- a/src/tearing.jl +++ b/src/tearing.jl @@ -49,3 +49,71 @@ function contract_variables(graph::BipartiteGraph, var_eq_matching::Matching, return newgraph end + +""" + $TYPEDSIGNATURES + +Find the equations (source vertices of `graph`) which are not matched to a variable present +in `vars_scc` according to the matching `var_eq_matching`. `varfilter` filters out +variables to exclude from this process. +""" +function free_equations(graph::BipartiteGraph, vars_scc::Vector{Vector{Int}}, + var_eq_matching::Matching, varfilter::F) where {F} + ne = nsrcs(graph) + seen_eqs = falses(ne) + for vars in vars_scc, var in vars + + varfilter(var) || continue + ieq = var_eq_matching[var] + if ieq isa Int + seen_eqs[ieq] = true + end + end + findall(!, seen_eqs) +end + +const MatchingT{T} = Matching{T, Vector{Union{T, Int}}} +const MatchedVarT = Union{Unassigned, SelectedState} +const VarEqMatchingT = MatchingT{MatchedVarT} + +""" + $TYPEDEF + +A struct containing the results of tearing. + +# Fields + +$TYPEDFIELDS +""" +struct TearingResult + """ + The variable-equation matching. Differential variables are matched to `SelectedState`. + The derivative of a differential variable is matched to the corresponding differential + equation. Solved variables are matched to the equation they are solved from. Algebraic + variables are matched to `unassigned`. + """ + var_eq_matching::VarEqMatchingT + """ + The variable-equation matching prior to tearing. This is the maximal matching used to + compute `var_sccs` (see below). For generating the torn system, `var_eq_matching` is + the source of truth. This should only be used to identify algebraic equations in each + SCC. + """ + full_var_eq_matching::VarEqMatchingT + """ + The partitioning of variables into strongly connected components (SCCs). The SCCs are + sorted in dependency order, so each SCC depends on variables in previous SCCs. + """ + var_sccs::Vector{Vector{Int}} +end + +""" + $TYPEDEF + +Supertype for all tearing algorithms. A tearing algorithm takes as input the +`SystemStructure` along with any other necessary arguments. + +The output of a tearing algorithm must be a `TearingResult` and a `NamedTuple` of +any additional data computed in the process that may be useful for further processing. +""" +abstract type TearingAlgorithm end From 956b3a2baf09e05aa0a2a261812414b0ae2e1ea3 Mon Sep 17 00:00:00 2001 From: Aayush Sabharwal Date: Tue, 2 Dec 2025 16:36:31 +0530 Subject: [PATCH 11/15] feat: add `trivial_tearing!` --- Project.toml | 2 + src/StateSelection.jl | 1 + src/tearing.jl | 106 ++++++++++++++++++++++++++++++++++++++++++ 3 files changed, 109 insertions(+) diff --git a/Project.toml b/Project.toml index ffe51ad..19929f3 100644 --- a/Project.toml +++ b/Project.toml @@ -9,6 +9,7 @@ DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae" FindFirstFunctions = "64ca27bc-2ba2-4a57-88aa-44e436879224" Graphs = "86223c79-3864-5bf0-83f7-82e725a168b6" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" +OrderedCollections = "bac558e1-5e72-5ebc-8fee-abe8a469f55d" Setfield = "efcf1570-3423-57d1-acb7-fd33fddbac46" SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" UnPack = "3a884ed6-31ef-47d7-9d2a-63182c4928ed" @@ -25,6 +26,7 @@ DocStringExtensions = "0.9.3" FindFirstFunctions = "1.2.0" Graphs = "1.10.0" LinearAlgebra = "1.11.0" +OrderedCollections = "1" Setfield = "1.1.1" SparseArrays = "1.11.0" UnPack = "1.0.2" diff --git a/src/StateSelection.jl b/src/StateSelection.jl index 3977b7c..dad40e2 100644 --- a/src/StateSelection.jl +++ b/src/StateSelection.jl @@ -5,6 +5,7 @@ using Setfield: @set!, @set using UnPack: @unpack using Graphs import SparseArrays +import OrderedCollections: OrderedSet # Graph Types using BipartiteGraphs diff --git a/src/tearing.jl b/src/tearing.jl index 02c1f15..bfe5bfb 100644 --- a/src/tearing.jl +++ b/src/tearing.jl @@ -50,6 +50,112 @@ function contract_variables(graph::BipartiteGraph, var_eq_matching::Matching, return newgraph end +""" + $(TYPEDSIGNATURES) + +Preemptively identify observed equations in the system and tear them. This happens before +any simplification. The equations torn by this process are ones that are already given in +an explicit form in the system and where the LHS is not present in any other equation of +the system except for other such preempitvely torn equations. +""" +function trivial_tearing!(ts::TransformationState) + # equations that can be trivially torn an observed equations + trivial_idxs = OrderedSet{Int}() + # variables that have been matched to trivially torn equations + matched_vars = OrderedSet{Int}() + + complete!(ts.structure) + var_to_diff = ts.structure.var_to_diff + graph = ts.structure.graph + candidates = collect(possibly_explicit_equations(ts)) + # TODO: Use DiCMOBiGraph here and topsort the equations. It'll remove the `while true`. + while true + # track whether we added an equation to the trivial list this iteration + added_equation = false + for (i, vari) in candidates + # don't check already torn equations + i in trivial_idxs && continue + + # if a variable was the LHS of two trivial observed equations, we wouldn't have + # included it in the list. Error if somehow it made it through. + @assert !(vari in matched_vars) + # don't tear differential/shift equations (or differentiated/shifted variables) + var_to_diff[vari] === nothing || continue + invview(var_to_diff)[vari] === nothing || continue + # get the equations that the candidate matched variable is present in, except + # those equations which have already been torn as observed + eqidxs = setdiff(𝑑neighbors(graph, vari), trivial_idxs) + # it should only be present in this equation + length(eqidxs) == 1 || continue + eqi = only(eqidxs) + @assert eqi == i + + # for every variable present in this equation, make sure it isn't _only_ + # present in trivial equations + isvalid = true + for v in 𝑠neighbors(graph, eqi) + v == vari && continue + v in matched_vars && continue + # `> 1` and not `0` because one entry will be this equation (`eqi`) + isvalid &= count(!in(trivial_idxs), 𝑑neighbors(graph, v)) > 1 + isvalid || break + end + isvalid || continue + + added_equation = true + push!(trivial_idxs, eqi) + push!(matched_vars, vari) + end + + # if we didn't add an equation this iteration, we won't add one next iteration + added_equation || break + end + + deleteat!(var_to_diff.primal_to_diff, matched_vars) + deleteat!(var_to_diff.diff_to_primal, matched_vars) + deleteat!(ts.structure.eq_to_diff.primal_to_diff, trivial_idxs) + deleteat!(ts.structure.eq_to_diff.diff_to_primal, trivial_idxs) + delete_srcs!(ts.structure.graph, trivial_idxs; rm_verts = true) + delete_dsts!(ts.structure.graph, matched_vars; rm_verts = true) + if ts.structure.solvable_graph !== nothing + delete_srcs!(ts.structure.solvable_graph, trivial_idxs; rm_verts = true) + delete_dsts!(ts.structure.solvable_graph, matched_vars; rm_verts = true) + end + trivial_tearing_postprocess!(ts, trivial_idxs, matched_vars) + return ts +end + +""" + $TYPEDSIGNATURES + +Return an iterable of tuples. The first element of each tuple is the index of an equation +index in `state` which has a single variable (present in `get_fullvars(state)`) on the LHS. +The second element of each tuple is the index of the variable on the LHS. + +These are considered candidates for [`trivial_tearing!`](@ref). Some equations may +intentionally be filtered out from this list, such as if the variable on the LHS should be +considered "irreducible" (not to be torn) or redundant equations which reduce to `0 ~ 0`. +""" +function possibly_explicit_equations(state::TransformationState) + error("This function must be implemented to run `trivial_tearing!`") +end + +""" + $TYPEDSIGNATURES + +Postprocessing function after running [`trivial_tearing!`](@ref). Update `state` given that +`torn_eqs` have been preemptively torn. The order of `torn_eqs` is important, as it +determines a topolgical ordering of the torn equations. `torn_vars` similarly identifies +the torn variables. The ordering of `torn_vars` corresponds to that of `torn_eqs`. + +Prior to calling this function, the minimal required fields of `state.structure` will have +been updated appropriately (torn elements removed). At minimum, this function should update +`state` such that [`get_fullvars`](@ref) returns the appropriate subset of variables. +""" +function trivial_tearing_postprocess!(state::TransformationState, torn_eqs::OrderedSet{Int}, torn_vars::OrderedSet{Int}) + error("This function must be implemented to run `trivial_tearing!`") +end + """ $TYPEDSIGNATURES From fdc9a1f995b7e41c45dd6406ace47d56be8e0811 Mon Sep 17 00:00:00 2001 From: Aayush Sabharwal Date: Tue, 2 Dec 2025 16:36:59 +0530 Subject: [PATCH 12/15] feat: add `ModiaTearing` algorithm --- src/StateSelection.jl | 2 +- src/modia_tearing.jl | 14 ++++++++++++++ 2 files changed, 15 insertions(+), 1 deletion(-) diff --git a/src/StateSelection.jl b/src/StateSelection.jl index dad40e2..63396b3 100644 --- a/src/StateSelection.jl +++ b/src/StateSelection.jl @@ -22,8 +22,8 @@ include("interface.jl") # Structural transformation passes include("singularity_removal.jl") include("pantelides.jl") -include("modia_tearing.jl") include("tearing.jl") +include("modia_tearing.jl") include("partial_state_selection.jl") # Utilities diff --git a/src/modia_tearing.jl b/src/modia_tearing.jl index a3c346f..aa6b17f 100644 --- a/src/modia_tearing.jl +++ b/src/modia_tearing.jl @@ -62,6 +62,20 @@ function tear_graph_block_modia!(var_eq_matching, ict, solvable_graph, eqs, vars return nothing end +@kwdef struct ModiaTearing{F, F2, F3} <: TearingAlgorithm + isder::F = nothing + varfilter::F2 = (_ -> true) + eqfilter::F3 = (_ -> true) +end + +function (alg::ModiaTearing)(structure::SystemStructure) + result = StateSelection.tear_graph_modia(structure, alg.isder; + varfilter = alg.varfilter, + eqfilter = alg.eqfilter) + var_eq_matching, full_var_eq_matching, var_sccs = result + return TearingResult(var_eq_matching, full_var_eq_matching, var_sccs), (;) +end + function tear_graph_modia(structure::SystemStructure, isder::F = nothing, ::Type{U} = Unassigned; varfilter::F2 = v -> true, From 6a633ef26b86b69b14d931ddec2e41244502406e Mon Sep 17 00:00:00 2001 From: Aayush Sabharwal Date: Tue, 2 Dec 2025 16:38:18 +0530 Subject: [PATCH 13/15] feat: make dummy derivatives implementation independent of MTK --- src/partial_state_selection.jl | 130 +++++++++++++++++++++++++++------ 1 file changed, 106 insertions(+), 24 deletions(-) diff --git a/src/partial_state_selection.jl b/src/partial_state_selection.jl index 5fc420d..0ac5764 100644 --- a/src/partial_state_selection.jl +++ b/src/partial_state_selection.jl @@ -176,60 +176,113 @@ end function dummy_derivative_graph!(state::TransformationState, jac = nothing; state_priority = nothing, log = Val(false), kwargs...) + state.structure.solvable_graph === nothing && find_solvables!(state; kwargs...) complete!(state.structure) var_eq_matching = complete(pantelides!(state)) dummy_derivative_graph!(state.structure, var_eq_matching, jac, state_priority, log) end +struct DummyDerivativeSummary + var_sccs::Vector{Vector{Int}} + state_priority::Vector{Vector{Float64}} +end + +""" + $TYPEDSIGNATURES + +Perform the dummy derivatives algorithm. + +# Arguments + +- `jac` is a function taking a list of equation and variable indices, and returning the + jacobian for the same if it has all integer entries. Otherwise, the function should + return `nothing`. +- `state_priority` is a function taking the index of a variable and returning its + priority. Higher priority variables are more likely to be chosen as states. +""" function dummy_derivative_graph!( structure::SystemStructure, var_eq_matching, jac = nothing, - state_priority = nothing, ::Val{log} = Val(false)) where {log} + state_priority = nothing, ::Val{log} = Val(false); + tearing_alg::TearingAlgorithm = DummyDerivativeTearing(), kwargs...) where {log} @unpack eq_to_diff, var_to_diff, graph = structure diff_to_eq = invview(eq_to_diff) diff_to_var = invview(var_to_diff) invgraph = invview(graph) + extended_sp = let state_priority = state_priority, var_to_diff = var_to_diff, + diff_to_var = diff_to_var + + var -> begin + min_p = max_p = 0.0 + while var_to_diff[var] !== nothing + var = var_to_diff[var] + end + while true + p = state_priority(var) + max_p = max(max_p, p) + min_p = min(min_p, p) + (var = diff_to_var[var]) === nothing && break + end + min_p < 0 ? min_p : max_p + end + end var_sccs = find_var_sccs(graph, var_eq_matching) + var_perm = Int[] + var_dummy_scc = Vector{Int}[] + var_state_priority = Vector{Float64}[] eqcolor = falses(nsrcs(graph)) dummy_derivatives = Int[] col_order = Int[] + neqs = nsrcs(graph) nvars = ndsts(graph) eqs = Int[] + vars = Int[] next_eq_idxs = Int[] next_var_idxs = Int[] new_eqs = Int[] new_vars = Int[] eqs_set = BitSet() - for vars in var_sccs + for vars′ in var_sccs empty!(eqs) - for var in vars + empty!(vars) + for var in vars′ eq = var_eq_matching[var] eq isa Int || continue - diff_to_eq[eq] === nothing && continue - push!(eqs, eq) + diff_to_eq[eq] === nothing || push!(eqs, eq) + if var_to_diff[var] !== nothing + error("Invalid SCC") + end + (diff_to_var[var] !== nothing && is_present(structure, var)) && push!(vars, var) end isempty(eqs) && continue - rank_matching = Matching(nvars) + rank_matching = Matching(max(nvars, neqs)) isfirst = true - if jac === nothing - J = nothing - else - _J = jac(eqs, vars) - # only accept small integers to avoid overflow - is_all_small_int = all(_J) do x′ - x = unwrap(x′) - x isa Number || return false - isinteger(x) && typemin(Int8) <= x <= typemax(Int8) + J = nothing + if jac !== nothing + J = jac(eqs, vars) + end + if J !== nothing + is_all_small_int = true + for i in eachindex(J) + x = J[i] + is_all_small_int = isinteger(x) && typemin(Int8) <= Int(x) <= typemax(Int8) + is_all_small_int || break end - J = is_all_small_int ? Int.(unwrap.(_J)) : nothing + J = is_all_small_int ? Int.(J) : nothing end while true nrows = length(eqs) iszero(nrows) && break if state_priority !== nothing && isfirst - sort!(vars, by = state_priority) + sp = extended_sp.(vars) + resize!(var_perm, length(sp)) + sortperm!(var_perm, sp) + permute!(vars, var_perm) + permute!(sp, var_perm) + push!(var_dummy_scc, copy(vars)) + push!(var_state_priority, sp) end # TODO: making the algorithm more robust # 1. If the Jacobian is a integer matrix, use Bareiss to check @@ -243,7 +296,7 @@ function dummy_derivative_graph!( if !isfirst J = J[next_eq_idxs, next_var_idxs] end - N = ModelingToolkit.nullspace(J; col_order) # modifies col_order + N = bareiss.nullspace(J; col_order) # modifies col_order rank = length(col_order) - size(N, 2) for i in 1:rank push!(dummy_derivatives, vars[col_order[i]]) @@ -291,6 +344,8 @@ function dummy_derivative_graph!( for (i, var) in enumerate(vars) ∫var = diff_to_var[var] ∫var === nothing && continue + ∫∫var = diff_to_var[∫var] + ∫∫var === nothing && continue if J !== nothing push!(next_var_idxs, i) end @@ -307,12 +362,9 @@ function dummy_derivative_graph!( @warn "The number of dummy derivatives ($n_dummys) does not match the number of differentiated equations ($n_diff_eqs)." end - ret = tearing_with_dummy_derivatives(structure, BitSet(dummy_derivatives)) - if log - ret - else - ret[1] - end + tearing_result, extra = tearing_alg(structure, BitSet(dummy_derivatives)) + extra = (; extra..., ddsummary = DummyDerivativeSummary(var_dummy_scc, var_state_priority)) + return tearing_result, extra end function is_present(structure, v)::Bool @@ -363,3 +415,33 @@ function tearing_with_dummy_derivatives(structure, dummy_derivatives) end return var_eq_matching, full_var_eq_matching, var_sccs, can_eliminate end + +struct DummyDerivativeTearing <: TearingAlgorithm end + +function (::DummyDerivativeTearing)(structure::SystemStructure, dummy_derivatives::Union{BitSet, Tuple{}} = ()) + @unpack var_to_diff = structure + # We can eliminate variables that are not selected (differential + # variables). Selected unknowns are differentiated variables that are not + # dummy derivatives. + can_eliminate = falses(length(var_to_diff)) + for (v, dv) in enumerate(var_to_diff) + dv = var_to_diff[v] + if dv === nothing || !is_some_diff(structure, dummy_derivatives, dv) + can_eliminate[v] = true + end + end + modia_tearing = ModiaTearing(; + isder = Base.Fix1(isdiffed, (structure, dummy_derivatives)), + varfilter = Base.Fix1(getindex, can_eliminate) + ) + tearing_result, _ = modia_tearing(structure) + + for v in 𝑑vertices(structure.graph) + is_present(structure, v) || continue + dv = var_to_diff[v] + (dv === nothing || !is_some_diff(structure, dummy_derivatives, dv)) && continue + tearing_result.var_eq_matching[v] = SelectedState() + end + + return tearing_result, (; can_eliminate) +end From 150ea7e5361799651f80350cea8db02c0f0df48b Mon Sep 17 00:00:00 2001 From: Aayush Sabharwal Date: Tue, 2 Dec 2025 16:38:28 +0530 Subject: [PATCH 14/15] chore: export some public API --- src/StateSelection.jl | 3 +++ 1 file changed, 3 insertions(+) diff --git a/src/StateSelection.jl b/src/StateSelection.jl index 63396b3..6f7bcc0 100644 --- a/src/StateSelection.jl +++ b/src/StateSelection.jl @@ -30,4 +30,7 @@ include("partial_state_selection.jl") include("debug.jl") include("utils.jl") +export DiffGraph, bareiss, CLIL, find_eq_solvables!, SelectedState +export TearingAlgorithm, TearingResult, ModiaTearing, DummyDerivativeTearing + end From 9dd2017bd5b245ca08090d38c9d768668296d27d Mon Sep 17 00:00:00 2001 From: Aayush Sabharwal Date: Tue, 2 Dec 2025 16:38:50 +0530 Subject: [PATCH 15/15] build: bump BipartiteGraphs compat --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 19929f3..2bd55b6 100644 --- a/Project.toml +++ b/Project.toml @@ -21,7 +21,7 @@ DeepDiffs = "ab62b9b5-e342-54a8-a765-a90f495de1a6" StateSelectionDeepDiffsExt = "DeepDiffs" [compat] -BipartiteGraphs = "0.1.1" +BipartiteGraphs = "0.1.2" DocStringExtensions = "0.9.3" FindFirstFunctions = "1.2.0" Graphs = "1.10.0"