Format

Author: vyudu

src/structural_transformation/StructuralTransformations.jl

@@ -11,7 +11,8 @@ using SymbolicUtils: maketerm, iscall
 
 using ModelingToolkit
 using ModelingToolkit: ODESystem, AbstractSystem, var_from_nested_derivative, Differential,
-                       unknowns, equations, vars, Symbolic, diff2term_with_unit, shift2term_with_unit, value,
+                       unknowns, equations, vars, Symbolic, diff2term_with_unit,
+                       shift2term_with_unit, value,
                        operation, arguments, Sym, Term, simplify, symbolic_linear_solve,
                        isdiffeq, isdifferential, isirreducible,
                        empty_substitutions, get_substitutions,
@@ -22,7 +23,8 @@ using ModelingToolkit: ODESystem, AbstractSystem, var_from_nested_derivative, Di
                        get_postprocess_fbody, vars!,
                        IncrementalCycleTracker, add_edge_checked!, topological_sort,
                        invalidate_cache!, Substitutions, get_or_construct_tearing_state,
-                       filter_kwargs, lower_varname_with_unit, lower_shift_varname_with_unit, setio, SparseMatrixCLIL,
+                       filter_kwargs, lower_varname_with_unit,
+                       lower_shift_varname_with_unit, setio, SparseMatrixCLIL,
                        get_fullvars, has_equations, observed,
                        Schedule, schedule
 

src/structural_transformation/symbolics_tearing.jl

@@ -248,7 +248,8 @@ called dummy derivatives.
 State selection is done. All non-differentiated variables are algebraic 
 variables, and all variables that appear differentiated are differential variables.
 """
-function substitute_derivatives_algevars!(ts::TearingState, neweqs, var_eq_matching, dummy_sub; iv = nothing, D = nothing)
+function substitute_derivatives_algevars!(
+        ts::TearingState, neweqs, var_eq_matching, dummy_sub; iv = nothing, D = nothing)
     @unpack fullvars, sys, structure = ts
     @unpack solvable_graph, var_to_diff, eq_to_diff, graph = structure
     diff_to_var = invview(var_to_diff)
@@ -288,7 +289,7 @@ end
 #= 
 There are three cases where we want to generate new variables to convert
 the system into first order (semi-implicit) ODEs.
-    
+
 1. To first order:
 Whenever higher order differentiated variable like `D(D(D(x)))` appears,
 we introduce new variables `x_t`, `x_tt`, and `x_ttt` and new equations
@@ -364,7 +365,8 @@ Effects on the system structure:
 - solvable_graph:
 - var_eq_matching: match D(x) to the added identity equation D(x) ~ x_t
 """
-function generate_derivative_variables!(ts::TearingState, neweqs, var_eq_matching; mm = nothing, iv = nothing, D = nothing)
+function generate_derivative_variables!(
+        ts::TearingState, neweqs, var_eq_matching; mm = nothing, iv = nothing, D = nothing)
     @unpack fullvars, sys, structure = ts
     @unpack solvable_graph, var_to_diff, eq_to_diff, graph = structure
     eq_var_matching = invview(var_eq_matching)
@@ -395,7 +397,7 @@ function generate_derivative_variables!(ts::TearingState, neweqs, var_eq_matchin
         dx = fullvars[dv]
         order, lv = var_order(dv, diff_to_var)
         x_t = is_discrete ? lower_shift_varname_with_unit(fullvars[dv], iv) :
-                          lower_varname_with_unit(fullvars[lv], iv, order)
+              lower_varname_with_unit(fullvars[lv], iv, order)
 
         # Add `x_t` to the graph
         v_t = add_dd_variable!(structure, fullvars, x_t, dv)
@@ -405,7 +407,7 @@ function generate_derivative_variables!(ts::TearingState, neweqs, var_eq_matchin
         # Update matching
         push!(var_eq_matching, unassigned)
         var_eq_matching[dv] = unassigned
-        eq_var_matching[dummy_eq] = dv 
+        eq_var_matching[dummy_eq] = dv
     end
 end
 
@@ -428,7 +430,7 @@ function find_duplicate_dd(dv, solvable_graph, diff_to_var, linear_eqs, mm)
             return eq, v_t
         end
     end
-    return nothing 
+    return nothing
 end
 
 """
@@ -492,8 +494,9 @@ Order the new equations and variables such that the differential equations
 and variables come first. Return the new equations, the solved equations,
 the new orderings, and the number of solved variables and equations.
 """
-function generate_system_equations!(state::TearingState, neweqs, var_eq_matching; simplify = false, iv = nothing, D = nothing)
-    @unpack fullvars, sys, structure = state 
+function generate_system_equations!(state::TearingState, neweqs, var_eq_matching;
+        simplify = false, iv = nothing, D = nothing)
+    @unpack fullvars, sys, structure = state
     @unpack solvable_graph, var_to_diff, eq_to_diff, graph = structure
     eq_var_matching = invview(var_eq_matching)
     diff_to_var = invview(var_to_diff)
@@ -502,11 +505,12 @@ function generate_system_equations!(state::TearingState, neweqs, var_eq_matching
     if is_only_discrete(structure)
         for (i, v) in enumerate(fullvars)
             op = operation(v)
-            op isa Shift && (op.steps < 0) && begin
-                lowered = lower_shift_varname_with_unit(v, iv)
-                total_sub[v] = lowered
-                fullvars[i] = lowered
-            end
+            op isa Shift && (op.steps < 0) &&
+                begin
+                    lowered = lower_shift_varname_with_unit(v, iv)
+                    total_sub[v] = lowered
+                    fullvars[i] = lowered
+                end
         end
     end
 
@@ -581,10 +585,11 @@ function generate_system_equations!(state::TearingState, neweqs, var_eq_matching
     end
     solved_vars_set = BitSet(solved_vars)
     var_ordering = [diff_vars;
-                   setdiff!(setdiff(1:ndsts(graph), diff_vars_set),
-                       solved_vars_set)]
+                    setdiff!(setdiff(1:ndsts(graph), diff_vars_set),
+                        solved_vars_set)]
 
-    return neweqs, solved_eqs, eq_ordering, var_ordering, length(solved_vars), length(solved_vars_set)
+    return neweqs, solved_eqs, eq_ordering, var_ordering, length(solved_vars),
+    length(solved_vars_set)
 end
 
 """
@@ -648,7 +653,8 @@ Eliminate the solved variables and equations from the graph and permute the
 graph's vertices to account for the new variable/equation ordering.
 """
 # TODO: BLT sorting
-function reorder_vars!(state::TearingState, var_eq_matching, eq_ordering, var_ordering, nsolved_eq, nsolved_var)
+function reorder_vars!(state::TearingState, var_eq_matching, eq_ordering,
+        var_ordering, nsolved_eq, nsolved_var)
     @unpack solvable_graph, var_to_diff, eq_to_diff, graph = state.structure
 
     eqsperm = zeros(Int, nsrcs(graph))
@@ -692,7 +698,8 @@ end
 """
 Update the system equations, unknowns, and observables after simplification.
 """
-function update_simplified_system!(state::TearingState, neweqs, solved_eqs, dummy_sub, var_eq_matching, extra_unknowns; 
+function update_simplified_system!(
+        state::TearingState, neweqs, solved_eqs, dummy_sub, var_eq_matching, extra_unknowns;
         cse_hack = true, array_hack = true)
     @unpack solvable_graph, var_to_diff, eq_to_diff, graph = state.structure
     diff_to_var = invview(var_to_diff)
@@ -732,7 +739,6 @@ function update_simplified_system!(state::TearingState, neweqs, solved_eqs, dumm
     sys = schedule(sys)
 end
 
-    
 """
 Give the order of the variable indexed by dv.
 """
@@ -790,12 +796,14 @@ function tearing_reassemble(state::TearingState, var_eq_matching,
 
     generate_derivative_variables!(state, neweqs, var_eq_matching; mm, iv, D)
 
-    neweqs, solved_eqs, eq_ordering, var_ordering, nelim_eq, nelim_var = 
-        generate_system_equations!(state, neweqs, var_eq_matching; simplify, iv, D)
+    neweqs, solved_eqs, eq_ordering, var_ordering, nelim_eq, nelim_var = generate_system_equations!(
+        state, neweqs, var_eq_matching; simplify, iv, D)
 
-    state = reorder_vars!(state, var_eq_matching, eq_ordering, var_ordering, nelim_eq, nelim_var)
+    state = reorder_vars!(
+        state, var_eq_matching, eq_ordering, var_ordering, nelim_eq, nelim_var)
 
-    sys = update_simplified_system!(state, neweqs, solved_eqs, dummy_sub, var_eq_matching, extra_unknowns; cse_hack, array_hack)
+    sys = update_simplified_system!(state, neweqs, solved_eqs, dummy_sub, var_eq_matching,
+        extra_unknowns; cse_hack, array_hack)
 
     @set! state.sys = sys
     @set! sys.tearing_state = state

src/structural_transformation/utils.jl

@@ -477,15 +477,17 @@ function shift2term(var)
     backshift = is_lowered ? op.steps + ModelingToolkit.getshift(arg) : op.steps
 
     num = join(Char(0x2080 + d) for d in reverse!(digits(-backshift))) # subscripted number, e.g. ₁
-    ds = join([Char(0x209c), Char(0x208b), num]) 
+    ds = join([Char(0x209c), Char(0x208b), num])
     # Char(0x209c) = ₜ
     # Char(0x208b) = ₋ (subscripted minus)
 
     O = is_lowered ? ModelingToolkit.getunshifted(arg) : arg
     oldop = operation(O)
-    newname = backshift != 0 ? Symbol(string(nameof(oldop)), ds) : Symbol(string(nameof(oldop)))
+    newname = backshift != 0 ? Symbol(string(nameof(oldop)), ds) :
+              Symbol(string(nameof(oldop)))
 
-    newvar = maketerm(typeof(O), Symbolics.rename(oldop, newname), Symbolics.children(O), Symbolics.metadata(O))
+    newvar = maketerm(typeof(O), Symbolics.rename(oldop, newname),
+        Symbolics.children(O), Symbolics.metadata(O))
     newvar = setmetadata(newvar, Symbolics.VariableSource, (:variables, newname))
     newvar = setmetadata(newvar, ModelingToolkit.VariableUnshifted, O)
     newvar = setmetadata(newvar, ModelingToolkit.VariableShift, backshift)

src/systems/discrete_system/discrete_system.jl

@@ -269,15 +269,15 @@ function shift_u0map_forward(sys::DiscreteSystem, u0map, defs)
     for k in collect(keys(u0map))
         v = u0map[k]
         if !((op = operation(k)) isa Shift)
-            isnothing(getunshifted(k)) && error("Initial conditions must be for the past state of the unknowns. Instead of providing the condition for $k, provide the condition for $(Shift(iv, -1)(k)).")
-            
+            isnothing(getunshifted(k)) &&
+                error("Initial conditions must be for the past state of the unknowns. Instead of providing the condition for $k, provide the condition for $(Shift(iv, -1)(k)).")
+
             updated[Shift(iv, 1)(k)] = v
         elseif op.steps > 0
             error("Initial conditions must be for the past state of the unknowns. Instead of providing the condition for $k, provide the condition for $(Shift(iv, -1)(only(arguments(k)))).")
         else
             updated[Shift(iv, op.steps + 1)(only(arguments(k)))] = v
         end
-
     end
     for var in unknowns(sys)
         op = operation(var)

src/systems/systemstructure.jl

@@ -473,9 +473,9 @@ function shift_discrete_system(ts::TearingState)
     end
     iv = get_iv(sys)
 
-    discmap = Dict(k => StructuralTransformations.simplify_shifts(Shift(iv, 1)(k))    
-                   for k in discvars 
-                   if any(isequal(k), fullvars) && !isa(operation(k), Union{Sample, Hold})) 
+    discmap = Dict(k => StructuralTransformations.simplify_shifts(Shift(iv, 1)(k))
+    for k in discvars
+    if any(isequal(k), fullvars) && !isa(operation(k), Union{Sample, Hold}))
 
     for i in eachindex(fullvars)
         fullvars[i] = StructuralTransformations.simplify_shifts(fast_substitute(

src/variables.jl

@@ -138,7 +138,7 @@ function default_toterm(x)
     if iscall(x) && (op = operation(x)) isa Operator
         if !(op isa Differential)
             if op isa Shift && op.steps < 0
-                return shift2term(x) 
+                return shift2term(x)
             end
             x = normalize_to_differential(op)(arguments(x)...)
         end

test/discrete_system.jl

@@ -257,7 +257,6 @@ k = ShiftIndex(t)
 @named sys = DiscreteSystem([x ~ x^2 + y^2, y ~ x(k - 1) + y(k - 1)], t)
 @test_throws ["algebraic equations", "not yet supported"] structural_simplify(sys)
 
-
 @testset "Passing `nothing` to `u0`" begin
     @variables x(t) = 1
     k = ShiftIndex()
@@ -273,11 +272,11 @@ end
     prob = DiscreteProblem(de, [], (0, 10))
     @test prob[x] == 2.0
     @test prob[x(k - 1)] == 1.0
-    
+
     # must provide initial conditions for history
     @test_throws ErrorException DiscreteProblem(de, [x => 2.0], (0, 10))
-    @test_throws ErrorException DiscreteProblem(de, [x(k+1) => 2.], (0, 10))
-    
+    @test_throws ErrorException DiscreteProblem(de, [x(k + 1) => 2.0], (0, 10))
+
     # initial values only affect _that timestep_, not the entire history
     prob = DiscreteProblem(de, [x(k - 1) => 2.0], (0, 10))
     @test prob[x] == 3.0
@@ -286,34 +285,35 @@ end
     @test prob[xₜ₋₁] == 2.0
 
     # Test initial assignment with lowered variable
-    prob = DiscreteProblem(de, [xₜ₋₁(k-1) => 4.0], (0, 10))
-    @test prob[x(k-1)] == prob[xₜ₋₁] == 1.0
-    @test prob[x] == 5.
+    prob = DiscreteProblem(de, [xₜ₋₁(k - 1) => 4.0], (0, 10))
+    @test prob[x(k - 1)] == prob[xₜ₋₁] == 1.0
+    @test prob[x] == 5.0
 
     # Test missing initial throws error
     @variables x(t)
-    @mtkbuild de = DiscreteSystem([x ~ x(k-1) + x(k-2)*x(k-3)], t)
-    @test_throws ErrorException prob = DiscreteProblem(de, [x(k-3) => 2.], (0, 10))
-    @test_throws ErrorException prob = DiscreteProblem(de, [x(k-3) => 2., x(k-1) => 3.], (0, 10))
+    @mtkbuild de = DiscreteSystem([x ~ x(k - 1) + x(k - 2) * x(k - 3)], t)
+    @test_throws ErrorException prob=DiscreteProblem(de, [x(k - 3) => 2.0], (0, 10))
+    @test_throws ErrorException prob=DiscreteProblem(
+        de, [x(k - 3) => 2.0, x(k - 1) => 3.0], (0, 10))
 
     # Test non-assigned initials are given default value
-    @variables x(t) = 2.
-    @mtkbuild de = DiscreteSystem([x ~ x(k-1) + x(k-2)*x(k-3)], t)
-    prob = DiscreteProblem(de, [x(k-3) => 12.], (0, 10))
+    @variables x(t) = 2.0
+    @mtkbuild de = DiscreteSystem([x ~ x(k - 1) + x(k - 2) * x(k - 3)], t)
+    prob = DiscreteProblem(de, [x(k - 3) => 12.0], (0, 10))
     @test prob[x] == 26.0
-    @test prob[x(k-1)] == 2.0
-    @test prob[x(k-2)] == 2.0
+    @test prob[x(k - 1)] == 2.0
+    @test prob[x(k - 2)] == 2.0
 
     # Elaborate test
     @variables xₜ₋₂(t) zₜ₋₁(t) z(t)
-    eqs = [x ~ x(k-1) + z(k-2), 
-           z ~ x(k-2) * x(k-3) - z(k-1)^2]
+    eqs = [x ~ x(k - 1) + z(k - 2),
+        z ~ x(k - 2) * x(k - 3) - z(k - 1)^2]
     @mtkbuild de = DiscreteSystem(eqs, t)
-    u0 = [x(k-1) => 3, 
-          xₜ₋₂(k-1) => 4, 
-          x(k-2) => 1, 
-          z(k-1) => 5, 
-          zₜ₋₁(k-1) => 12]
+    u0 = [x(k - 1) => 3,
+        xₜ₋₂(k - 1) => 4,
+        x(k - 2) => 1,
+        z(k - 1) => 5,
+        zₜ₋₁(k - 1) => 12]
     prob = DiscreteProblem(de, u0, (0, 10))
     @test prob[x] == 15
     @test prob[z] == -21