From 2b1bf2f382b65ac112b92d2ff853914941fc76c8 Mon Sep 17 00:00:00 2001 From: Yingbo Ma Date: Wed, 3 Feb 2021 21:43:22 -0500 Subject: [PATCH 1/5] Add `check` kw in `solve_for` --- src/solve.jl | 21 ++++++++++++--------- 1 file changed, 12 insertions(+), 9 deletions(-) diff --git a/src/solve.jl b/src/solve.jl index 206f63638d..a8faceba52 100644 --- a/src/solve.jl +++ b/src/solve.jl @@ -56,34 +56,37 @@ end # return the coefficient matrix `A` and a # vector of constants (possibly symbolic) `b` such that # A \ b will solve the equations for the vars -function A_b(eqs::AbstractArray, vars::AbstractArray) +function A_b(eqs::AbstractArray, vars::AbstractArray, check) exprs = rhss(eqs) .- lhss(eqs) - for ex in exprs - @assert islinear(ex, vars) + if check + for ex in exprs + @assert islinear(ex, vars) + end end A = jacobian(exprs, vars) b = A * vars - exprs A, b end -function A_b(eq, var) +function A_b(eq, var, check) ex = eq.rhs - eq.lhs - @assert islinear(ex, [var]) + check && @assert islinear(ex, [var]) a = expand_derivatives(Differential(var)(ex)) b = a * var - ex a, b end """ - solve_for(eqs::Vector, vars::Vector) + solve_for(eqs::Vector, vars::Vector; simplify=true, check=true) Solve the vector of equations `eqs` for a set of variables `vars`. Assumes `length(eqs) == length(vars)` -Currently only works if all equations are linear. +Currently only works if all equations are linear. `check` if the expr is linear +w.r.t `vars`. """ -function solve_for(eqs, vars; simplify=true) - A, b = A_b(eqs, vars) +function solve_for(eqs, vars; simplify=true, check=true) + A, b = A_b(eqs, vars, check) #TODO: we need to make sure that `solve_for(eqs, vars)` contains no `vars` _solve(A, b, simplify) end From 3c4ea7d379bf4436e99947708167968b632698c3 Mon Sep 17 00:00:00 2001 From: Yingbo Ma Date: Thu, 4 Feb 2021 01:37:21 -0500 Subject: [PATCH 2/5] Add some convenient interfaces --- src/systems/systemstructure.jl | 49 +++++++++++++++++++++++++++++++--- 1 file changed, 45 insertions(+), 4 deletions(-) diff --git a/src/systems/systemstructure.jl b/src/systems/systemstructure.jl index 2f6e95d723..e38e7d7ea3 100644 --- a/src/systems/systemstructure.jl +++ b/src/systems/systemstructure.jl @@ -1,3 +1,11 @@ +module SystemStructures + +using ..ModelingToolkit +using ..ModelingToolkit: isdiffeq, var_from_nested_derivative, vars! +using SymbolicUtils: arguments +using ..BipartiteGraphs +using UnPack +using Setfield using SparseArrays #= @@ -27,10 +35,18 @@ for v in š‘£vertices(graph); active_š‘£vertices[v] || continue end =# +export SystemStructure, initialize_system_structure +export diffvars_range, dervars_range, algvars_range +export isdiffvars, isdervars, isalgvars +export DIFFERENTIAL_VARIABLE, ALGEBRAIC_VARIABLE, DERIVATIVE_VARIABLE +export DIFFERENTIAL_EQUATION, ALGEBRAIC_EQUATION +export vartype, eqtype + struct SystemStructure dxvar_offset::Int fullvars::Vector # [xvar; dxvars; algvars] varassoc::Vector{Int} + algeqs::BitVector graph::BipartiteGraph{Int} solvable_graph::BipartiteGraph{Int} assign::Vector{Int} @@ -39,12 +55,33 @@ struct SystemStructure partitions::Vector{NTuple{4, Vector{Int}}} end +diffvars_range(s::SystemStructure) = 1:s.dxvar_offset +dervars_range(s::SystemStructure) = s.dxvar_offset+1:2s.dxvar_offset +algvars_range(s::SystemStructure) = 2s.dxvar_offset+1:length(s.fullvars) + +isdiffvar(s::SystemStructure, var::Integer) = var in diffvars_range(s) +isdervar(s::SystemStructure, var::Integer) = var in dervars_range(s) +isalgvar(s::SystemStructure, var::Integer) = var in algvars_range(s) + +@enum VariableType DIFFERENTIAL_VARIABLE ALGEBRAIC_VARIABLE DERIVATIVE_VARIABLE + +function vartype(s::SystemStructure, var::Integer)::VariableType + isdiffvar(s, var) ? DIFFERENTIAL_VARIABLE : + isdervar(s, var) ? DERIVATIVE_VARIABLE : + isalgvar(s, var) ? ALGEBRAIC_VARIABLE : error("Variable $var out of bounds") +end + +@enum EquationType DIFFERENTIAL_EQUATION ALGEBRAIC_EQUATION + +eqtype(s::SystemStructure, eq::Integer)::EquationType = s.algeqs[eq] ? ALGEBRAIC_EQUATION : DIFFERENTIAL_EQUATION + function initialize_system_structure(sys) - sys, dxvar_offset, fullvars, varassoc, graph, solvable_graph = init_graph(flatten(sys)) + sys, dxvar_offset, fullvars, varassoc, algeqs, graph, solvable_graph = init_graph(sys) @set sys.structure = SystemStructure( dxvar_offset, fullvars, varassoc, + algeqs, graph, solvable_graph, Int[], @@ -74,8 +111,10 @@ end function collect_variables(sys) dxvars = [] eqs = equations(sys) + algeqs = falses(length(eqs)) for (i, eq) in enumerate(eqs) if isdiffeq(eq) + algeqs[i] = true lhs = eq.lhs # Make sure that the LHS is a first order derivative of a var. @assert !(arguments(lhs)[1] isa Differential) "The equation $eq is not first order" @@ -86,11 +125,11 @@ function collect_variables(sys) xvars = (first āˆ˜ var_from_nested_derivative).(dxvars) algvars = setdiff(states(sys), xvars) - return xvars, dxvars, algvars + return xvars, dxvars, algvars, algeqs end function init_graph(sys) - xvars, dxvars, algvars = collect_variables(sys) + xvars, dxvars, algvars, algeqs = collect_variables(sys) dxvar_offset = length(xvars) algvar_offset = 2dxvar_offset @@ -119,5 +158,7 @@ function init_graph(sys) end varassoc = Int[(1:dxvar_offset) .+ dxvar_offset; zeros(Int, length(fullvars) - dxvar_offset)] # variable association list - sys, dxvar_offset, fullvars, varassoc, graph, solvable_graph + sys, dxvar_offset, fullvars, varassoc, algeqs, graph, solvable_graph end + +end # module From 9dcf55f8c3b5719f2af71e72d71473492923df91 Mon Sep 17 00:00:00 2001 From: Yingbo Ma Date: Thu, 4 Feb 2021 01:37:55 -0500 Subject: [PATCH 3/5] Modularize BipartiteGraphs and SystemStructures --- Project.toml | 2 ++ src/ModelingToolkit.jl | 6 +++++- src/bipartite_graph.jl | 13 ++++++++++++- 3 files changed, 19 insertions(+), 2 deletions(-) diff --git a/Project.toml b/Project.toml index 778427b0f1..42afb9a634 100644 --- a/Project.toml +++ b/Project.toml @@ -21,6 +21,7 @@ LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" MacroTools = "1914dd2f-81c6-5fcd-8719-6d5c9610ff09" NaNMath = "77ba4419-2d1f-58cd-9bb1-8ffee604a2e3" RecursiveArrayTools = "731186ca-8d62-57ce-b412-fbd966d074cd" +Reexport = "189a3867-3050-52da-a836-e630ba90ab69" Requires = "ae029012-a4dd-5104-9daa-d747884805df" RuntimeGeneratedFunctions = "7e49a35a-f44a-4d26-94aa-eba1b4ca6b47" SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f" @@ -49,6 +50,7 @@ LightGraphs = "1.3" MacroTools = "0.5" NaNMath = "0.3" RecursiveArrayTools = "2.3" +Reexport = "1" Requires = "1.0" RuntimeGeneratedFunctions = "0.4, 0.5" SafeTestsets = "0.0.1" diff --git a/src/ModelingToolkit.jl b/src/ModelingToolkit.jl index c5f96c85b9..fdef0b884a 100644 --- a/src/ModelingToolkit.jl +++ b/src/ModelingToolkit.jl @@ -5,6 +5,7 @@ using StaticArrays, LinearAlgebra, SparseArrays, LabelledArrays using Latexify, Unitful, ArrayInterface using MacroTools using UnPack: @unpack +using Setfield using DiffEqJump using DataStructures using SpecialFunctions, NaNMath @@ -202,6 +203,7 @@ Get the set of parameters variables for the given system. function parameters end include("bipartite_graph.jl") +using .BipartiteGraphs include("variables.jl") include("context_dsl.jl") @@ -216,7 +218,6 @@ include("domains.jl") include("register_function.jl") include("systems/abstractsystem.jl") -include("systems/systemstructure.jl") include("systems/diffeqs/odesystem.jl") include("systems/diffeqs/sdesystem.jl") @@ -239,6 +240,9 @@ include("systems/pde/pdesystem.jl") include("systems/reaction/reactionsystem.jl") include("systems/dependency_graphs.jl") +include("systems/systemstructure.jl") +using .SystemStructures + include("systems/reduction.jl") include("latexify_recipes.jl") diff --git a/src/bipartite_graph.jl b/src/bipartite_graph.jl index c8cea617a7..5d5ea83d19 100644 --- a/src/bipartite_graph.jl +++ b/src/bipartite_graph.jl @@ -1,6 +1,15 @@ +module BipartiteGraphs + +export BipartiteEdge, BipartiteGraph + +export š‘ vertices, š‘‘vertices, has_š‘ vertex, has_š‘‘vertex, š‘ neighbors, š‘‘neighbors, + š‘ edges, š‘‘edges, nsrcs, ndsts, SRC, DST + +using DocStringExtensions +using Reexport using UnPack using SparseArrays -using LightGraphs +@reexport using LightGraphs using Setfield ### @@ -229,3 +238,5 @@ function LightGraphs.incidence_matrix(g::BipartiteGraph, val=true) end S = sparse(I, J, val, nsrcs(g), ndsts(g)) end + +end # module From a78fbf19e229a7271535b962fde60eef3317692a Mon Sep 17 00:00:00 2001 From: Yingbo Ma Date: Thu, 4 Feb 2021 01:39:41 -0500 Subject: [PATCH 4/5] Fix a typo --- src/systems/systemstructure.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/systems/systemstructure.jl b/src/systems/systemstructure.jl index e38e7d7ea3..16babf1c0b 100644 --- a/src/systems/systemstructure.jl +++ b/src/systems/systemstructure.jl @@ -1,7 +1,7 @@ module SystemStructures using ..ModelingToolkit -using ..ModelingToolkit: isdiffeq, var_from_nested_derivative, vars! +using ..ModelingToolkit: isdiffeq, var_from_nested_derivative, vars!, flatten using SymbolicUtils: arguments using ..BipartiteGraphs using UnPack @@ -76,7 +76,7 @@ end eqtype(s::SystemStructure, eq::Integer)::EquationType = s.algeqs[eq] ? ALGEBRAIC_EQUATION : DIFFERENTIAL_EQUATION function initialize_system_structure(sys) - sys, dxvar_offset, fullvars, varassoc, algeqs, graph, solvable_graph = init_graph(sys) + sys, dxvar_offset, fullvars, varassoc, algeqs, graph, solvable_graph = init_graph(flatten(sys)) @set sys.structure = SystemStructure( dxvar_offset, fullvars, From a85f70200ccc869aa90797e8c407e1fb9bdcae93 Mon Sep 17 00:00:00 2001 From: Yingbo Ma Date: Thu, 4 Feb 2021 10:59:27 -0500 Subject: [PATCH 5/5] Add methods isdiffeq and isalgeq on a SystemStructure --- src/systems/systemstructure.jl | 12 +++++++----- 1 file changed, 7 insertions(+), 5 deletions(-) diff --git a/src/systems/systemstructure.jl b/src/systems/systemstructure.jl index 16babf1c0b..a475458c80 100644 --- a/src/systems/systemstructure.jl +++ b/src/systems/systemstructure.jl @@ -1,7 +1,7 @@ module SystemStructures using ..ModelingToolkit -using ..ModelingToolkit: isdiffeq, var_from_nested_derivative, vars!, flatten +import ..ModelingToolkit: isdiffeq, var_from_nested_derivative, vars!, flatten using SymbolicUtils: arguments using ..BipartiteGraphs using UnPack @@ -37,7 +37,7 @@ end export SystemStructure, initialize_system_structure export diffvars_range, dervars_range, algvars_range -export isdiffvars, isdervars, isalgvars +export isdiffvar, isdervar, isalgvar, isdiffeq, isalgeq export DIFFERENTIAL_VARIABLE, ALGEBRAIC_VARIABLE, DERIVATIVE_VARIABLE export DIFFERENTIAL_EQUATION, ALGEBRAIC_EQUATION export vartype, eqtype @@ -73,7 +73,9 @@ end @enum EquationType DIFFERENTIAL_EQUATION ALGEBRAIC_EQUATION -eqtype(s::SystemStructure, eq::Integer)::EquationType = s.algeqs[eq] ? ALGEBRAIC_EQUATION : DIFFERENTIAL_EQUATION +isalgeq(s::SystemStructure, eq::Integer) = s.algeqs[eq] +isdiffeq(s::SystemStructure, eq::Integer) = !isalgeq(s, eq) +eqtype(s::SystemStructure, eq::Integer)::EquationType = isalgeq(s, eq) ? ALGEBRAIC_EQUATION : DIFFERENTIAL_EQUATION function initialize_system_structure(sys) sys, dxvar_offset, fullvars, varassoc, algeqs, graph, solvable_graph = init_graph(flatten(sys)) @@ -111,10 +113,10 @@ end function collect_variables(sys) dxvars = [] eqs = equations(sys) - algeqs = falses(length(eqs)) + algeqs = trues(length(eqs)) for (i, eq) in enumerate(eqs) if isdiffeq(eq) - algeqs[i] = true + algeqs[i] = false lhs = eq.lhs # Make sure that the LHS is a first order derivative of a var. @assert !(arguments(lhs)[1] isa Differential) "The equation $eq is not first order"