Import MTK code

LICENSE

@@ -1,6 +1,6 @@
 MIT License
 
-Copyright (c) 2024 JuliaHub, Inc. and other contributors
+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

Project.toml

@@ -1,11 +1,33 @@
 name = "StateSelection"
 uuid = "64909d44-ed92-46a8-bbd9-f047dfbdc84b"
+version = "0.1.0-DEV"
 authors = ["JuliaHub", "Inc. and other contributors"]
-version = "1.0.0-DEV"
+
+[deps]
+DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
+FindFirstFunctions = "64ca27bc-2ba2-4a57-88aa-44e436879224"
+Graphs = "86223c79-3864-5bf0-83f7-82e725a168b6"
+LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+Setfield = "efcf1570-3423-57d1-acb7-fd33fddbac46"
+SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
+UnPack = "3a884ed6-31ef-47d7-9d2a-63182c4928ed"
 
 [compat]
+DocStringExtensions = "0.9.3"
+FindFirstFunctions = "1.2.0"
+Graphs = "1.10.0"
+LinearAlgebra = "1.11.0"
+Setfield = "1.1.1"
+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"
 

README.md

@@ -1,5 +1,79 @@
 # StateSelection
 
-[![Build Status](https://github.com/Keno/StateSelection.jl/actions/workflows/CI.yml/badge.svg?branch=main)](https://github.com/Keno/StateSelection.jl/actions/workflows/CI.yml?query=branch%3Amain)
-[![Coverage](https://codecov.io/gh/Keno/StateSelection.jl/branch/main/graph/badge.svg)](https://codecov.io/gh/Keno/StateSelection.jl)
-[![Coverage](https://coveralls.io/repos/github/Keno/StateSelection.jl/badge.svg?branch=main)](https://coveralls.io/github/Keno/StateSelection.jl?branch=main)
+[![Build Status](https://github.com/JuliaComputing/StateSelection.jl/actions/workflows/CI.yml/badge.svg?branch=main)](https://github.com/JuliaComputing/StateSelection.jl/actions/workflows/CI.yml?query=branch%3Amain)
+[![Coverage](https://codecov.io/gh/JuliaComputing/StateSelection.jl/branch/main/graph/badge.svg)](https://codecov.io/gh/JuliaComputing/StateSelection.jl)
+[![Coverage](https://coveralls.io/repos/github/JuliaComputing/StateSelection.jl/badge.svg?branch=main)](https://coveralls.io/github/JuliaComputing/StateSelection.jl?branch=main)
+
+This package implements *structural* transformations suitable
+of optimizing systems of (non-linear, ordinary differential, differential algebraic) equations for faster and more stable integration using a numerical solver. It is intended to serve as a common algorithmic core to a variety of downstream modeling systems, including [MTK](https://github.com/SciML/ModelingToolkit.jl), [DAECompiler](https://github.com/CedarEDA/DAECompiler.jl) and JuliaSimCompiler.
+
+It is intended to be independent of the actual representation of the system of equations, instead operating on structural abstractions defined in this package. In particular, it computes only the transformations that *should* be the done to the system of equations, but the actual transformation of the system itself is deferred to the downstream user of this package.
+
+## Transformations
+
+This package implements *partial state selection* PSS. Unfortunately, the terminology is not entirely consistent in the literature, so before we proceed, we define PSS problem.
+
+### The state selection problem
+
+We will consider a *system* to be the dictum of a list of variables and equations together with
+
+1. The structural incidence matrix of equations and variables
+2. The numerical incidence matrix of the integer-linear constant-coefficient sub-system
+3. For each equation-variable pair an indication of whether the downstream compiler is capable of solving the equation for that variable.
+4. A graph of differential relations between the variables
+
+Then, the *PSS problem* is to find
+
+A. A (possibly empty) list of equations to differentiate
+B. The assignment of each variable (or their derivatives if such derivative occurs in a differentiated equation) to one of
+
+ i. A selected differential state
+ ii. An assignment to an equation (declaring the equation will be solved for this variable)
+ iii. An assignment to a linear-system of a list of equations
+ iv. An algebraic state
+
+Such that
+1. the dependency graph of variables (described further below) is acyclic
+2. The provided solvability constraints are obeyed.
+
+### Further remarks
+
+There are several equivalent ways of thinking of this assignment. Two common ones in the literature are
+1. an upper-triangular structural incidence matrix
+2. a graph of dependencies between variables
+
+To form the dependency graph of variables, let the vertices of the graph be the variables. For each variable assigned to an equation, add a directed edge from all other variables in the incidence-row of the said equation to the assigned variable.
+
+Solutions to the PSS problem are not unique and the chosen solution materially affects the ease and stability of the resulting numerical integration. In general, a numerical integrator may need to switch the set of chosen states dynamically at runtime (known as dynamic state selection). Note that while this package has various hooks to control the selected states, and can compute the tearing sets, it does not by itself implement dynamic state selection.
+
+# Detailed description of the transformation steps
+
+A. Structural Singularity Removal
+
+This is a pre-pass that runs before pantelides. The primary objective of the pantelides algorithm is to ensure that the
+jacobian of the fundamental ODE system is non-singular at runtime. However, because pantelides is a structural algorithm,
+it can accidentally fail to fully reduce the index for systems which have full structural rank, but are numerically singular.
+One common source of such systems are conservation laws commonly used in hierachical modeling tools.
+
+However, fortunately, in such systems, the numerical singularity is static and apparent in the constant linear coefficients of the integer-linear subsystem (ILS) of the whole system.
+
+As such, this pre-pass applies a change of basis to the ILS to make the numerical rank-deficiency structurally apparent, allowing pantelides to properly differentiate the resulting equations.
+
+B. Patenlides algorithm (DAE only)
+
+The Pantelides algorithm is used for reduxing the differentiation index of the system to index 1 or 0, making it suitable for integration with a numerical integrator (such integrators are generally not capable of integrating higher-index DAE systems). This is accomplished by generating a list of equations to differentiate and adding these to the system.
+
+C. Dummy Derivatives (DAE only)
+
+The pantelides algorithm produces systems that are over-determined in their differential relations (the sum of the number of differential relations and algebraic relations needs to match the number of variables). As such, some of these differential relations need to be removed in order to make the system solvable. However, the choice of which relations to remove can radically affect the numerical properties of the system and thus the ease and stability of the resulting integration.
+
+D. Tearing
+
+Tearing computes a (directed) dependency graph of equations (or dually variables). If this graph is fully connected and has no cycles, all output variables are uniquely determined given all input variables, so the system will have no algebraic states. In general however, this dependency graph may have cycles (referred to as *algebraic loops*). Such cycles are broken numerically by choosing one or more variables in the cycle as algebraic states and (at runtime) wrapping a non-linear solver around these variables. Again,
+the choice of algebraic states greatly affects the ease and stability of the resulting solve or integration.
+
+
+## History
+
+The code in this package can be traced back to
+JuliaComputing/StructuralTransformations.jl, which provided a port of the structural transformation core from [Modia](https://github.com/ModiaSim/Modia.jl) to MTK's data structures. This package was subsequently integrated into MTK, with a number of intervening rewrites, simplifications and cleanups (both in code and understanding), before be-ing re-extracted to become this package.

ext/StateSelectionDeepDiffsExt.jl

@@ -0,0 +1,187 @@
+
+using DeepDiffs, ModelingToolkit
+using StateSelection.BipartiteGraphs: Label,
+ BipartiteAdjacencyList, unassigned,
+ HighlightInt
+using StateSelection: SystemStructure,
+ MatchedSystemStructure,
+ SystemStructurePrintMatrix
+
+"""
+A utility struct for displaying the difference between two HighlightInts.
+
+# Example
+```julia
+using StateSelection, DeepDiffs
+
+old_i = HighlightInt(1, :default, true)
+new_i = HighlightInt(2, :default, false)
+diff = HighlightIntDiff(new_i, old_i)
+
+show(diff)
+```
+"""
+struct HighlightIntDiff
+ new::HighlightInt
+ old::HighlightInt
+end
+
+function Base.show(io::IO, d::HighlightIntDiff)
+ p_color = d.new.highlight
+ (d.new.match && !d.old.match) && (p_color = :light_green)
+ (!d.new.match && d.old.match) && (p_color = :light_red)
+
+ (d.new.match || d.old.match) && printstyled(io, "(", color = p_color)
+ if d.new.i != d.old.i
+ Base.show(io, HighlightInt(d.old.i, :light_red, d.old.match))
+ print(io, " ")
+ Base.show(io, HighlightInt(d.new.i, :light_green, d.new.match))
+ else
+ Base.show(io, HighlightInt(d.new.i, d.new.highlight, false))
+ end
+ (d.new.match || d.old.match) && printstyled(io, ")", color = p_color)
+end
+
+"""
+A utility struct for displaying the difference between two
+BipartiteAdjacencyList's.
+
+# Example
+```julia
+using ModelingToolkit, DeepDiffs
+
+old = BipartiteAdjacencyList(...)
+new = BipartiteAdjacencyList(...)
+diff = BipartiteAdjacencyListDiff(new, old)
+
+show(diff)
+```
+"""
+struct BipartiteAdjacencyListDiff
+ new::BipartiteAdjacencyList
+ old::BipartiteAdjacencyList
+end
+
+function Base.show(io::IO, l::BipartiteAdjacencyListDiff)
+ print(io,
+ LabelDiff(Label(l.new.match === true ? "∫ " : ""),
+ Label(l.old.match === true ? "∫ " : "")))
+ (l.new.match !== true && l.old.match !== true) && print(io, " ")
+
+ new_nonempty = isnothing(l.new.u) ? nothing : !isempty(l.new.u)
+ old_nonempty = isnothing(l.old.u) ? nothing : !isempty(l.old.u)
+ if new_nonempty === true && old_nonempty === true
+ if (!isempty(setdiff(l.new.highlight_u, l.new.u)) ||
+ !isempty(setdiff(l.old.highlight_u, l.old.u)))
+ throw(ArgumentError("The provided `highlight_u` must be a sub-graph of `u`."))
+ end
+
+ new_items = Dict(i => HighlightInt(i, :nothing, i === l.new.match) for i in l.new.u)
+ old_items = Dict(i => HighlightInt(i, :nothing, i === l.old.match) for i in l.old.u)
+
+ highlighted = union(map(intersect(l.new.u, l.old.u)) do i
+ HighlightIntDiff(new_items[i], old_items[i])
+ end,
+ map(setdiff(l.new.u, l.old.u)) do i
+ HighlightInt(new_items[i].i, :light_green,
+ new_items[i].match)
+ end,
+ map(setdiff(l.old.u, l.new.u)) do i
+ HighlightInt(old_items[i].i, :light_red,
+ old_items[i].match)
+ end)
+ print(IOContext(io, :typeinfo => typeof(highlighted)), highlighted)
+ elseif new_nonempty === true
+ printstyled(
+ io, map(l.new.u) do i
+ HighlightInt(i, :nothing, i === l.new.match)
+ end, color = :light_green)
+ elseif old_nonempty === true
+ printstyled(
+ io, map(l.old.u) do i
+ HighlightInt(i, :nothing, i === l.old.match)
+ end, color = :light_red)
+ elseif old_nonempty !== nothing || new_nonempty !== nothing
+ print(io,
+ LabelDiff(Label(new_nonempty === false ? "∅" : "", :light_black),
+ Label(old_nonempty === false ? "∅" : "", :light_black)))
+ else
+ printstyled(io, '⋅', color = :light_black)
+ end
+end
+
+"""
+A utility struct for displaying the difference between two Labels
+in git-style red/green highlighting.
+
+# Example
+```julia
+using ModelingToolkit, DeepDiffs
+
+old = Label("before")
+new = Label("after")
+diff = LabelDiff(new, old)
+
+show(diff)
+```
+"""
+struct LabelDiff
+ new::Label
+ old::Label
+end
+function Base.show(io::IO, l::LabelDiff)
+ if l.new != l.old
+ printstyled(io, l.old.s, color = :light_red)
+ length(l.new.s) != 0 && length(l.old.s) != 0 && print(io, " ")
+ printstyled(io, l.new.s, color = :light_green)
+ else
+ print(io, l.new)
+ end
+end
+
+"""
+A utility struct for displaying the difference between two
+(Matched)SystemStructure's in git-style red/green highlighting.
+
+# Example
+```julia
+using ModelingToolkit, DeepDiffs
+
+old = SystemStructurePrintMatrix(...)
+new = SystemStructurePrintMatrix(...)
+diff = SystemStructureDiffPrintMatrix(new, old)
+
+show(diff)
+```
+"""
+struct SystemStructureDiffPrintMatrix <:
+ AbstractMatrix{Union{LabelDiff, BipartiteAdjacencyListDiff}}
+ new::SystemStructurePrintMatrix
+ old::SystemStructurePrintMatrix
+end
+
+function Base.size(ssdpm::SystemStructureDiffPrintMatrix)
+ max.(Base.size(ssdpm.new), Base.size(ssdpm.old))
+end
+
+function Base.getindex(ssdpm::SystemStructureDiffPrintMatrix, i::Integer, j::Integer)
+ checkbounds(ssdpm, i, j)
+ if i > 1 && (j == 4 || j == 9)
+ old = new = BipartiteAdjacencyList(nothing, nothing, unassigned)
+ (i <= size(ssdpm.new, 1)) && (new = ssdpm.new[i, j])
+ (i <= size(ssdpm.old, 1)) && (old = ssdpm.old[i, j])
+ BipartiteAdjacencyListDiff(new, old)
+ else
+ old = new = Label("")
+ (i <= size(ssdpm.new, 1)) && (new = ssdpm.new[i, j])
+ (i <= size(ssdpm.old, 1)) && (old = ssdpm.old[i, j])
+ LabelDiff(new, old)
+ end
+end
+
+function DeepDiffs.deepdiff(old::Union{MatchedSystemStructure, SystemStructure},
+ new::Union{MatchedSystemStructure, SystemStructure})
+ new_sspm = SystemStructurePrintMatrix(new)
+ old_sspm = SystemStructurePrintMatrix(old)
+ Base.print_matrix(stdout, SystemStructureDiffPrintMatrix(new_sspm, old_sspm))
+end

src/StateSelection.jl

@@ -1,5 +1,34 @@
 module StateSelection
 
-# Write your package code here.
+using DocStringExtensions
+using Setfield: @set!, @set
+using UnPack: @unpack
+using Graphs
+
+# Graph Types
+function invview end
+function complete end
+include("graph/bipartite.jl")
+include("graph/diff.jl")
+using .BipartiteGraphs
+
+# Math library
+include("math/bareiss.jl")
+include("math/sparsematrixclil.jl")
+using .CLIL, .bareiss
+
+# Downstream interferace
+include("interface.jl")
+
+# Structural transformation passes
+include("singularity_removal.jl")
+include("pantelides.jl")
+include("modia_tearing.jl")
+include("tearing.jl")
+include("partial_state_selection.jl")
+
+# Utilities
+include("debug.jl")
+include("utils.jl")
 
 end