diff --git a/lib/ModelingToolkitBase/src/ModelingToolkitBase.jl b/lib/ModelingToolkitBase/src/ModelingToolkitBase.jl index 2ccd665935..69db941a71 100644 --- a/lib/ModelingToolkitBase/src/ModelingToolkitBase.jl +++ b/lib/ModelingToolkitBase/src/ModelingToolkitBase.jl @@ -10,6 +10,9 @@ using PrecompileTools, Reexport using JumpProcesses # ONLY here for the invalidations import REPL + using OffsetArrays: Origin + import BlockArrays: BlockArray, BlockedArray, Block, blocksize, blocksizes, blockpush!, + undef_blocks, blocks end import SymbolicUtils @@ -60,9 +63,6 @@ using Moshi.Data: @data using Reexport using RecursiveArrayTools import Graphs: SimpleDiGraph, add_edge!, incidence_matrix -import BlockArrays: BlockArray, BlockedArray, Block, blocksize, blocksizes, blockpush!, - undef_blocks, blocks -using OffsetArrays: Origin import CommonSolve import EnumX import ReadOnlyDicts: ReadOnlyDict diff --git a/lib/ModelingToolkitBase/src/discretedomain.jl b/lib/ModelingToolkitBase/src/discretedomain.jl index 6531596340..c005ffff51 100644 --- a/lib/ModelingToolkitBase/src/discretedomain.jl +++ b/lib/ModelingToolkitBase/src/discretedomain.jl @@ -23,7 +23,7 @@ struct Shift <: Operator """Fixed Shift""" t::Union{Nothing, SymbolicT} steps::Int - Shift(t, steps = 1) = new(value(t), steps) + Shift(t, steps = 1) = new(unwrap(t), steps) end Shift(steps::Int) = new(nothing, steps) normalize_to_differential(s::Shift) = Differential(s.t)^s.steps diff --git a/lib/ModelingToolkitBase/src/precompile.jl b/lib/ModelingToolkitBase/src/precompile.jl index 26cc55a574..9369304b69 100644 --- a/lib/ModelingToolkitBase/src/precompile.jl +++ b/lib/ModelingToolkitBase/src/precompile.jl @@ -12,6 +12,7 @@ PrecompileTools.@compile_workload begin x ^ 5 6 ^ x x - y + x * x * q[1] -y 2y z = 2 @@ -81,11 +82,22 @@ PrecompileTools.@compile_workload begin q[1] q'q using ModelingToolkitBase - @variables x(ModelingToolkitBase.t_nounits) y(ModelingToolkitBase.t_nounits) - isequal(ModelingToolkitBase.D_nounits.x, ModelingToolkitBase.t_nounits) + @parameters g + @variables x(ModelingToolkitBase.t_nounits) + @variables y(ModelingToolkitBase.t_nounits) [state_priority = 10] + @variables λ(ModelingToolkitBase.t_nounits) + eqs = [ + ModelingToolkitBase.D_nounits(ModelingToolkitBase.D_nounits(x)) ~ λ * x + ModelingToolkitBase.D_nounits(ModelingToolkitBase.D_nounits(y)) ~ λ * y - g + x^2 + y^2 ~ 1 + ] + dvs = Num[x, y, λ] + ps = Num[g] ics = Dict{SymbolicT, SymbolicT}() - ics[x] = 2.3 - sys = System([ModelingToolkitBase.D_nounits(x) ~ x * y, y ~ 2x], ModelingToolkitBase.t_nounits, [x, y], Num[]; initial_conditions = ics, guesses = ics, name = :sys) + ics[y] = -1.0 + ics[ModelingToolkitBase.D_nounits(x)] = 0.5 + isequal(ModelingToolkitBase.D_nounits.x, ModelingToolkitBase.t_nounits) + sys = System(eqs, ModelingToolkitBase.t_nounits, dvs, ps; initial_conditions = ics, guesses = ics, name = :sys) complete(sys) @static if @isdefined(ModelingToolkit) TearingState(sys) diff --git a/lib/ModelingToolkitBase/src/problems/jumpproblem.jl b/lib/ModelingToolkitBase/src/problems/jumpproblem.jl index 03b8fd2aa6..ad29b16677 100644 --- a/lib/ModelingToolkitBase/src/problems/jumpproblem.jl +++ b/lib/ModelingToolkitBase/src/problems/jumpproblem.jl @@ -172,12 +172,11 @@ end ##### MTK dispatches for Symbolic jumps ##### eqtype_supports_collect_vars(j::MassActionJump) = true -function collect_vars!(unknowns::OrderedSet{SymbolicT}, parameters::OrderedSet{SymbolicT}, j::MassActionJump, iv::Union{SymbolicT, Nothing}; depth = 0, - op = Differential) - collect_vars!(unknowns, parameters, j.scaled_rates, iv; depth, op) +function collect_vars!(unknowns::OrderedSet{SymbolicT}, parameters::OrderedSet{SymbolicT}, j::MassActionJump, iv::Union{SymbolicT, Nothing}, ::Type{op} = Differential; depth = 0) where {op} + collect_vars!(unknowns, parameters, j.scaled_rates, iv, op; depth) for field in (j.reactant_stoch, j.net_stoch) for el in field - collect_vars!(unknowns, parameters, el, iv; depth, op) + collect_vars!(unknowns, parameters, el, iv, op; depth) end end return nothing @@ -185,10 +184,10 @@ end eqtype_supports_collect_vars(j::Union{ConstantRateJump, VariableRateJump}) = true function collect_vars!(unknowns::OrderedSet{SymbolicT}, parameters::OrderedSet{SymbolicT}, j::Union{ConstantRateJump, VariableRateJump}, - iv::Union{SymbolicT, Nothing}; depth = 0, op = Differential) - collect_vars!(unknowns, parameters, j.rate, iv; depth, op) + iv::Union{SymbolicT, Nothing}, ::Type{op} = Differential; depth = 0) where {op} + collect_vars!(unknowns, parameters, j.rate, iv, op; depth) for eq in j.affect! - (eq isa Equation) && collect_vars!(unknowns, parameters, eq, iv; depth, op) + (eq isa Equation) && collect_vars!(unknowns, parameters, eq, iv, op; depth) end return nothing end diff --git a/lib/ModelingToolkitBase/src/systems/abstractsystem.jl b/lib/ModelingToolkitBase/src/systems/abstractsystem.jl index 880108ac58..7112da2c57 100644 --- a/lib/ModelingToolkitBase/src/systems/abstractsystem.jl +++ b/lib/ModelingToolkitBase/src/systems/abstractsystem.jl @@ -1460,6 +1460,14 @@ function unknowns_toplevel(sys::AbstractSystem) return get_unknowns(sys) end +function __no_initial_params_pred(x::SymbolicT) + arr, _ = split_indexed_var(x) + Moshi.Match.@match arr begin + BSImpl.Term(; f) && if f isa Initial end => false + _ => true + end +end + """ $(TYPEDSIGNATURES) @@ -1482,11 +1490,7 @@ function parameters(sys::AbstractSystem; initial_parameters = false) end result = collect(result) if !initial_parameters && !is_initializesystem(sys) - filter!(result) do sym - return !(isoperator(sym, Initial) || - iscall(sym) && operation(sym) === getindex && - isoperator(arguments(sym)[1], Initial)) - end + filter!(__no_initial_params_pred, result) end return result end @@ -1699,7 +1703,7 @@ Recursively substitute `dict` into `expr`. Use `Symbolics.simplify` on the expre if `simplify == true`. """ function substitute_and_simplify(expr, dict::AbstractDict, simplify::Bool) - expr = Symbolics.fixpoint_sub(expr, dict; operator = Union{ModelingToolkitBase.Initial, Pre}) + expr = Symbolics.fixpoint_sub(expr, dict, Union{Initial, Pre}) simplify ? Symbolics.simplify(expr) : expr end diff --git a/lib/ModelingToolkitBase/src/systems/callbacks.jl b/lib/ModelingToolkitBase/src/systems/callbacks.jl index 68a3e18256..3b7af2b5c3 100644 --- a/lib/ModelingToolkitBase/src/systems/callbacks.jl +++ b/lib/ModelingToolkitBase/src/systems/callbacks.jl @@ -89,7 +89,7 @@ function AffectSystem(affect::Vector{Equation}; discrete_parameters = SymbolicT[ if !haspre(eq) && !(isconst(eq.lhs) && isconst(eq.rhs)) @invokelatest warn_algebraic_equation(eq) end - collect_vars!(dvs, params, eq, iv; op = Pre) + collect_vars!(dvs, params, eq, iv, Pre) empty!(_varsbuf) SU.search_variables!(_varsbuf, eq; is_atomic = OperatorIsAtomic{Pre}()) filter!(x -> iscall(x) && operation(x) === Pre(), _varsbuf) @@ -125,7 +125,7 @@ function AffectSystem(affect::Vector{Equation}; discrete_parameters = SymbolicT[ # This `@invokelatest` should not be necessary, but it works around the inference bug # in https://github.com/JuliaLang/julia/issues/59943. Remove it at your own risk, the # bug took weeks to reduce to an MWE. - affectsys = @invokelatest mtkcompile(affectsys; fully_determined = nothing) + affectsys = (@invokelatest mtkcompile(affectsys; fully_determined = nothing))::System # get accessed parameters p from Pre(p) in the callback parameters accessed_params = Vector{SymbolicT}(filter(isparameter, map(unPre, collect(pre_params)))) union!(accessed_params, sys_params) @@ -133,12 +133,27 @@ function AffectSystem(affect::Vector{Equation}; discrete_parameters = SymbolicT[ # add scalarized unknowns to the map. _obs, _ = unhack_observed(observed(affectsys), equations(affectsys)) _dvs = vcat(unknowns(affectsys), map(eq -> eq.lhs, _obs)) - _dvs = reduce(vcat, map(safe_vec ∘ scalarize, _dvs), init = SymbolicT[]) - _discs = reduce(vcat, map(safe_vec ∘ scalarize, discretes); init = SymbolicT[]) + _dvs = __safe_scalarize_vars(_dvs) + _discs = __safe_scalarize_vars(discretes) setdiff!(_dvs, _discs) AffectSystem(affectsys, _dvs, accessed_params, discrete_parameters) end +function __safe_scalarize_vars(vars::Vector{SymbolicT}) + _vars = SymbolicT[] + for v in vars + sh = SU.shape(v)::SU.ShapeVecT + if isempty(sh) + push!(_vars, v) + continue + end + for i in SU.stable_eachindex(v) + push!(_vars, v[i]) + end + end + return _vars +end + safe_vec(@nospecialize(x)) = x isa SymbolicT ? [x] : vec(x::Array{SymbolicT}) system(a::AffectSystem) = a.system @@ -1043,13 +1058,13 @@ The `SymbolicDiscreteCallback`s in the returned vector are structs with two fiel See also `get_discrete_events`, which only returns the events of the top-level system. """ function discrete_events(sys::AbstractSystem) - obs = get_discrete_events(sys) + cbs = get_discrete_events(sys) systems = get_systems(sys) - cbs = [obs; - reduce(vcat, - (map(cb -> namespace_callback(cb, s), discrete_events(s)) for s in systems), - init = SymbolicDiscreteCallback[])] - cbs + cbs = copy(cbs) + for s in systems + append!(cbs, map(Base.Fix2(namespace_callback, s), discrete_events(s))) + end + return cbs end """ @@ -1100,15 +1115,13 @@ The `SymbolicContinuousCallback`s in the returned vector are structs with two fi See also `get_continuous_events`, which only returns the events of the top-level system. """ function continuous_events(sys::AbstractSystem) - obs = get_continuous_events(sys) - filter(!isempty, obs) - + cbs = get_continuous_events(sys) systems = get_systems(sys) - cbs = [obs; - reduce(vcat, - (map(o -> namespace_callback(o, s), continuous_events(s)) for s in systems), - init = SymbolicContinuousCallback[])] - filter(!isempty, cbs) + cbs = copy(cbs) + for s in systems + append!(cbs, map(Base.Fix2(namespace_callback, s), continuous_events(s))) + end + return cbs end """ diff --git a/lib/ModelingToolkitBase/src/systems/codegen_utils.jl b/lib/ModelingToolkitBase/src/systems/codegen_utils.jl index 320d0be769..92d9b2c884 100644 --- a/lib/ModelingToolkitBase/src/systems/codegen_utils.jl +++ b/lib/ModelingToolkitBase/src/systems/codegen_utils.jl @@ -42,7 +42,7 @@ reconstruct array variables if they are present scalarized in `args`. """ function array_variable_assignments(args...; argument_name = generated_argument_name) # map array symbolic to an identically sized array where each element is (buffer_idx, idx_in_buffer) - var_to_arridxs = Dict{BasicSymbolic, Array{Tuple{Int, Int}}}() + var_to_arridxs = Dict{SymbolicT, Vector{Tuple{Int, Int}}}() for (i, arg) in enumerate(args) # filter out non-arrays # any element of args which is not an array is assumed to not contain a @@ -55,13 +55,12 @@ function array_variable_assignments(args...; argument_name = generated_argument_ for (j, var) in enumerate(arg) var = unwrap(var) # filter out non-array-symbolics - iscall(var) || continue - operation(var) == getindex || continue - arrvar = arguments(var)[1] + arrvar, isarr = split_indexed_var(var) + isarr || continue # get and/or construct the buffer storing indexes idxbuffer = get!( - () -> map(Returns((0, 0)), eachindex(arrvar)), var_to_arridxs, arrvar) - Origin(first.(axes(arrvar))...)(idxbuffer)[arguments(var)[2:end]...] = (i, j) + () -> map(Returns((0, 0)), SU.stable_eachindex(arrvar)), var_to_arridxs, arrvar) + idxbuffer[SU.as_linear_idx(SU.shape(arrvar), get_stable_index(var))] = (i, j) end end diff --git a/lib/ModelingToolkitBase/src/systems/system.jl b/lib/ModelingToolkitBase/src/systems/system.jl index 83eff004fe..77a0f4ef37 100644 --- a/lib/ModelingToolkitBase/src/systems/system.jl +++ b/lib/ModelingToolkitBase/src/systems/system.jl @@ -381,12 +381,12 @@ function System(eqs::Vector{Equation}, iv, dvs, ps, brownians = SymbolicT[]; continuous_events = SymbolicContinuousCallback[], discrete_events = SymbolicDiscreteCallback[], connector_type = nothing, assertions = Dict{SymbolicT, String}(), metadata = MetadataT(), gui_metadata = nothing, - is_dde = nothing, tstops = [], inputs = OrderedSet{SymbolicT}(), + is_dde = nothing, @nospecialize(tstops = []), inputs = OrderedSet{SymbolicT}(), outputs = OrderedSet{SymbolicT}(), tearing_state = nothing, ignored_connections = nothing, parent = nothing, description = "", name = nothing, discover_from_metadata = true, initializesystem = nothing, is_initializesystem = false, is_discrete = false, - preface = [], checks = true, __legacy_defaults__ = nothing) + @nospecialize(preface = nothing), checks = true, __legacy_defaults__ = nothing) name === nothing && throw(NoNameError()) if __legacy_defaults__ !== nothing @@ -480,14 +480,17 @@ function System(eqs::Vector{Equation}, iv, dvs, ps, brownians = SymbolicT[]; filter!(!(Base.Fix1(===, COMMON_NOTHING) ∘ last), guesses) if iv === nothing - filter!(bindings) do kvp - k = kvp[1] - if k in all_dvs - initial_conditions[k] = kvp[2] - return false + filterer = let initial_conditions = initial_conditions, all_dvs = all_dvs + function _filterer(kvp) + k = kvp[1] + if k in all_dvs + initial_conditions[k] = kvp[2] + return false + end + return true end - return true end + filter!(filterer, bindings) end check_bindings(ps, bindings) @@ -683,19 +686,20 @@ Create a `System` with a single equation `eq`. System(eq::Equation, args...; kwargs...) = System([eq], args...; kwargs...) function gather_array_params(ps) - new_ps = OrderedSet() + new_ps = OrderedSet{SymbolicT}() for p in ps - if iscall(p) && operation(p) === getindex - par = arguments(p)[begin] - if symbolic_has_known_size(par) && all(par[i] in ps for i in eachindex(par)) - push!(new_ps, par) + arr, isarr = split_indexed_var(p) + sh = SU.shape(arr) + if isarr + if !(sh isa SU.Unknown) && all(in(ps) ∘ Base.Fix1(getindex, arr), SU.stable_eachindex(arr)) + push!(new_ps, arr) else push!(new_ps, p) end else - if symbolic_type(p) == ArraySymbolic() && symbolic_has_known_size(p) - for i in eachindex(p) - delete!(new_ps, p[i]) + if sh isa SU.ShapeVecT && !isempty(sh) + for i in SU.stable_eachindex(arr) + delete!(new_ps, arr[i]) end end push!(new_ps, p) diff --git a/lib/ModelingToolkitBase/src/systems/systems.jl b/lib/ModelingToolkitBase/src/systems/systems.jl index 1a2c007dba..69ec599a2e 100644 --- a/lib/ModelingToolkitBase/src/systems/systems.jl +++ b/lib/ModelingToolkitBase/src/systems/systems.jl @@ -10,14 +10,15 @@ function canonicalize_io(iovars, type::String) iobuffer = OrderedSet{SymbolicT}() arrsyms = AtomicArrayDict{OrderedSet{SymbolicT}}() for var in iovars - if Symbolics.isarraysymbolic(var) - if !symbolic_has_known_size(var) - throw(ArgumentError(""" - All $(type)s must have known shape. Found $var with unknown shape. - """)) + sh = SU.shape(var) + if SU.is_array_shape(sh) + if sh isa SU.ShapeVecT + union!(iobuffer, vec(collect(var)::Array{SymbolicT})::Vector{SymbolicT}) + continue end - union!(iobuffer, vec(collect(var)::Array{SymbolicT})::Vector{SymbolicT}) - continue + throw(ArgumentError(""" + All $(type)s must have known shape. Found $var with unknown shape. + """)) end arr, isarr = split_indexed_var(var) if isarr @@ -41,7 +42,7 @@ function canonicalize_io(iovars, type::String) or simply pass $k as an $type. """)) end - if type != "output" && !isequal(vec(collect(k))::Vector{SymbolicT}, collect(v)) + if type != "output" && !isequal(vec(collect(k)::Array{SymbolicT})::Vector{SymbolicT}, collect(v)) throw(ArgumentError(""" Elements of scalarized array variables must be in sorted order in $(type)s. \ Either pass all scalarized elements in sorted order as $(type)s \ @@ -469,16 +470,17 @@ function __num_isdiag_noise(mat) true end -function __get_num_diag_noise(mat) - map(axes(mat, 1)) do i +function __get_num_diag_noise(mat::Matrix{SymbolicT}) + result = fill(Symbolics.COMMON_ZERO, size(mat, 1)) + for i in axes(mat, 1) for j in axes(mat, 2) mij = mat[i, j] - if !_iszero(mij) - return mij - end + _iszero(mij) && continue + result[i] = mij + break end - 0 end + return result end """ diff --git a/lib/ModelingToolkitBase/src/utils.jl b/lib/ModelingToolkitBase/src/utils.jl index 4de348e2ed..833ef016ea 100644 --- a/lib/ModelingToolkitBase/src/utils.jl +++ b/lib/ModelingToolkitBase/src/utils.jl @@ -495,7 +495,7 @@ end """ Check if all the LHS are unique """ -function check_operator_variables(eqs, op::T) where {T} +function check_operator_variables(eqs, ::Type{op}) where {op} ops = Set{SymbolicT}() tmp = Set{SymbolicT}() for eq in eqs @@ -599,26 +599,26 @@ Search through equations and parameter dependencies of `sys`, where sys is at a recursively searches through all subsystems of `sys`, increasing the depth if it is not `-1`. A depth of `-1` indicates searching for variables with `GlobalScope`. """ -function collect_scoped_vars!(unknowns::OrderedSet{SymbolicT}, parameters::OrderedSet{SymbolicT}, sys::AbstractSystem, iv::Union{SymbolicT, Nothing}; depth = 1, op = Differential) +function collect_scoped_vars!(unknowns::OrderedSet{SymbolicT}, parameters::OrderedSet{SymbolicT}, sys::AbstractSystem, iv::Union{SymbolicT, Nothing}, ::Type{op} = Differential; depth = 1) where {op} if has_eqs(sys) for eq in equations(sys) eqtype_supports_collect_vars(eq) || continue if eq isa Equation symtype(eq.lhs) <: Number || continue end - collect_vars!(unknowns, parameters, eq, iv; depth, op) + collect_vars!(unknowns, parameters, eq, iv, op; depth) end end if has_jumps(sys) for eq in jumps(sys) eqtype_supports_collect_vars(eq) || continue - collect_vars!(unknowns, parameters, eq, iv; depth, op) + collect_vars!(unknowns, parameters, eq, iv, op; depth) end end if has_constraints(sys) for eq in constraints(sys) eqtype_supports_collect_vars(eq) || continue - collect_vars!(unknowns, parameters, eq, iv; depth, op) + collect_vars!(unknowns, parameters, eq, iv, op; depth) end end end @@ -710,7 +710,7 @@ can be checked using `check_scope_depth`. This function should return `nothing`. """ -function collect_vars!(unknowns::OrderedSet{SymbolicT}, parameters::OrderedSet{SymbolicT}, expr::SymbolicT, iv::Union{SymbolicT, Nothing}; depth = 0, op = Symbolics.Operator) +function collect_vars!(unknowns::OrderedSet{SymbolicT}, parameters::OrderedSet{SymbolicT}, expr::SymbolicT, iv::Union{SymbolicT, Nothing}, ::Type{op} = Symbolics.Operator; depth = 0) where {op} Moshi.Match.@match expr begin BSImpl.Const(;) => return BSImpl.Sym(;) => return collect_var!(unknowns, parameters, expr, iv; depth) @@ -731,9 +731,9 @@ function collect_vars!(unknowns::OrderedSet{SymbolicT}, parameters::OrderedSet{S return nothing end -function collect_vars!(unknowns::OrderedSet{SymbolicT}, parameters::OrderedSet{SymbolicT}, expr::AbstractArray, iv::Union{SymbolicT, Nothing}; depth = 0, op = Symbolics.Operator) +function collect_vars!(unknowns::OrderedSet{SymbolicT}, parameters::OrderedSet{SymbolicT}, expr::AbstractArray, iv::Union{SymbolicT, Nothing}, ::Type{op} = Symbolics.Operator; depth = 0) where {op} for var in expr - collect_vars!(unknowns, parameters, var, iv; depth, op) + collect_vars!(unknowns, parameters, var, iv, op; depth) end return nothing end @@ -749,27 +749,27 @@ eqtype_supports_collect_vars(eq::Equation) = true eqtype_supports_collect_vars(eq::Inequality) = true eqtype_supports_collect_vars(eq::Pair) = true -function collect_vars!(unknowns::OrderedSet{SymbolicT}, parameters::OrderedSet{SymbolicT}, eq::Union{Equation, Inequality}, iv::Union{SymbolicT, Nothing}; - depth = 0, op = Symbolics.Operator) - collect_vars!(unknowns, parameters, eq.lhs, iv; depth, op) - collect_vars!(unknowns, parameters, eq.rhs, iv; depth, op) +function collect_vars!(unknowns::OrderedSet{SymbolicT}, parameters::OrderedSet{SymbolicT}, eq::Union{Equation, Inequality}, iv::Union{SymbolicT, Nothing}, ::Type{op} = Symbolics.Operator; + depth = 0) where {op} + collect_vars!(unknowns, parameters, eq.lhs, iv, op; depth) + collect_vars!(unknowns, parameters, eq.rhs, iv, op; depth) return nothing end function collect_vars!( - unknowns::OrderedSet{SymbolicT}, parameters::OrderedSet{SymbolicT}, ex::Union{Num, Arr, CallAndWrap}, iv::Union{SymbolicT, Nothing}; depth = 0, op = Symbolics.Operator) - collect_vars!(unknowns, parameters, unwrap(ex), iv; depth, op) + unknowns::OrderedSet{SymbolicT}, parameters::OrderedSet{SymbolicT}, ex::Union{Num, Arr, CallAndWrap}, iv::Union{SymbolicT, Nothing}, ::Type{op} = Symbolics.Operator; depth = 0) where {op} + collect_vars!(unknowns, parameters, unwrap(ex), iv, op; depth) end function collect_vars!( - unknowns::OrderedSet{SymbolicT}, parameters::OrderedSet{SymbolicT}, p::Pair, iv::Union{SymbolicT, Nothing}; depth = 0, op = Symbolics.Operator) - collect_vars!(unknowns, parameters, p[1], iv; depth, op) - collect_vars!(unknowns, parameters, p[2], iv; depth, op) + unknowns::OrderedSet{SymbolicT}, parameters::OrderedSet{SymbolicT}, p::Pair, iv::Union{SymbolicT, Nothing}, ::Type{op} = Symbolics.Operator; depth = 0) where {op} + collect_vars!(unknowns, parameters, p[1], iv, op; depth) + collect_vars!(unknowns, parameters, p[2], iv, op; depth) return nothing end function collect_vars!( - unknowns::OrderedSet{SymbolicT}, parameters::OrderedSet{SymbolicT}, expr, iv::Union{SymbolicT, Nothing}; depth = 0, op = Symbolics.Operator) + unknowns::OrderedSet{SymbolicT}, parameters::OrderedSet{SymbolicT}, expr, iv::Union{SymbolicT, Nothing}, ::Type{op} = Symbolics.Operator; depth = 0) where {op} return nothing end @@ -806,7 +806,17 @@ function collect_var!(unknowns::OrderedSet{SymbolicT}, parameters::OrderedSet{Sy end # Add also any parameters that appear only as defaults in the var if hasdefault(var) && (def = getdefault(var)) !== missing - collect_vars!(unknowns, parameters, def, iv) + if def isa SymbolicT + collect_vars!(unknowns, parameters, def, iv) + elseif def isa Num + collect_vars!(unknowns, parameters, def, iv) + elseif def isa Arr + collect_vars!(unknowns, parameters, def, iv) + elseif def isa CallAndWrap + collect_vars!(unknowns, parameters, def, iv) + else + collect_vars!(unknowns, parameters, def, iv) + end end return nothing end diff --git a/lib/ModelingToolkitBase/src/variables.jl b/lib/ModelingToolkitBase/src/variables.jl index 16f811eb05..0bb7941e79 100644 --- a/lib/ModelingToolkitBase/src/variables.jl +++ b/lib/ModelingToolkitBase/src/variables.jl @@ -288,7 +288,7 @@ function shift2term(var::SymbolicT) newargs = copy(parent(args)) newargs[1] = shift2term(op(newargs[1])) unshifted_args = copy(newargs) - unshifted_args[1] = ModelingToolkitBase.getunshifted(newargs[1]) + unshifted_args[1] = ModelingToolkitBase.getunshifted(newargs[1])::SymbolicT unshifted = BSImpl.Term{VartypeT}(getindex, unshifted_args; type, shape, metadata) if metadata === nothing metadata = Base.ImmutableDict{DataType, Any}(VariableUnshifted, unshifted) @@ -302,7 +302,7 @@ function shift2term(var::SymbolicT) unshifted = ModelingToolkitBase.getunshifted(arg) is_lowered = unshifted !== nothing backshift = op.steps + ModelingToolkitBase.getshift(arg) - iszero(backshift) && return unshifted + iszero(backshift) && return unshifted::SymbolicT io = IOBuffer() O = (is_lowered ? unshifted : arg)::SymbolicT write(io, getname(O)) @@ -359,29 +359,24 @@ Simplify multiple shifts: Shift(t, k1)(Shift(t, k2)(x)) becomes Shift(t, k1+k2)( """ function simplify_shifts(var::SymbolicT) ModelingToolkitBase.hasshift(var) || return var - return SU.Rewriters.Postwalk(_simplify_shifts)(var) + return SU.Rewriters.Postwalk(_simplify_shifts)(var)::SymbolicT end +distribute_shift(eq::Equation) = distribute_shift(eq.lhs) ~ distribute_shift(eq.rhs) +distribute_shift(var::Union{Num, Arr}) = distribute_shift(unwrap(var)) """ Distribute a shift applied to a whole expression or equation. Shift(t, 1)(x + y) will become Shift(t, 1)(x) + Shift(t, 1)(y). Only shifts variables whose independent variable is the same t that appears in the Shift (i.e. constants, time-independent parameters, etc. do not get shifted). """ -function distribute_shift(var) - var = unwrap(var) - var isa Equation && return distribute_shift(var.lhs) ~ distribute_shift(var.rhs) - - ModelingToolkitBase.hasshift(var) || return var - shift = operation(var) - shift isa Shift || return var - - shift = operation(var) - expr = only(arguments(var)) - if expr isa Equation - return distribute_shift(shift(expr.lhs)) ~ distribute_shift(shift(expr.rhs)) +function distribute_shift(var::SymbolicT) + Moshi.Match.@match var begin + BSImpl.Term(; f, args) && if f isa Shift end => begin + shiftexpr = _distribute_shift(args[1], f) + return simplify_shifts(shiftexpr) + end + _ => return var end - shiftexpr = _distribute_shift(expr, shift) - return simplify_shifts(shiftexpr) end """ @@ -391,10 +386,10 @@ Whether `distribute_shift` should distribute shifts into the given operation. """ distribute_shift_into_operator(_) = true -function _distribute_shift(expr, shift) +function _distribute_shift(expr::SymbolicT, shift) if iscall(expr) op = operation(expr) - distribute_shift_into_operator(op) || return expr + distribute_shift_into_operator(op)::Bool || return expr args = arguments(expr) if ModelingToolkitBase.isvariable(expr) && operation(expr) !== getindex && diff --git a/lib/ModelingToolkitBase/test/analysis_points.jl b/lib/ModelingToolkitBase/test/analysis_points.jl index 08c3282a11..95d5afa616 100644 --- a/lib/ModelingToolkitBase/test/analysis_points.jl +++ b/lib/ModelingToolkitBase/test/analysis_points.jl @@ -6,6 +6,9 @@ using Test using ModelingToolkitBase: t_nounits as t, D_nounits as D, AnalysisPoint, AbstractSystem import ModelingToolkitBase as MTK import ControlSystemsBase as CS +using SciCompDSL +using ModelingToolkitStandardLibrary + using Symbolics: NAMESPACE_SEPARATOR @testset "AnalysisPoint is lowered to `connect`" begin diff --git a/lib/ModelingToolkitBase/test/symbolic_events.jl b/lib/ModelingToolkitBase/test/symbolic_events.jl index d65566a536..f13261ff54 100644 --- a/lib/ModelingToolkitBase/test/symbolic_events.jl +++ b/lib/ModelingToolkitBase/test/symbolic_events.jl @@ -646,7 +646,7 @@ end ps = @parameters k=k Θ=0.5 eqs = [D(x) ~ v, D(v) ~ -k * x + F] ev = [x ~ Θ] => [x ~ 1.0, v ~ 0.0] - System(eqs, t, sts, ps, continuous_events = [ev]; name) + System(eqs, t, sts, ps; continuous_events = [ev], name) end @named oscce = oscillator_ce() diff --git a/src/ModelingToolkit.jl b/src/ModelingToolkit.jl index 2818e3f0b8..3da86c06fb 100644 --- a/src/ModelingToolkit.jl +++ b/src/ModelingToolkit.jl @@ -8,6 +8,9 @@ using PrecompileTools, Reexport using Symbolics # ONLY here for the invalidations import REPL + using OffsetArrays: Origin + import BlockArrays: BlockArray, BlockedArray, Block, blocksize, blocksizes, blockpush!, + undef_blocks, blocks end import SymbolicUtils @@ -50,7 +53,6 @@ using Moshi.Data: @data import SCCNonlinearSolve using Reexport import Graphs: SimpleDiGraph, add_edge!, incidence_matrix -using OffsetArrays: Origin import CommonSolve using RuntimeGeneratedFunctions @@ -86,8 +88,6 @@ import PreallocationTools import PreallocationTools: DiffCache import FillArrays using BipartiteGraphs -import BlockArrays: BlockArray, BlockedArray, Block, blocksize, blocksizes, blockpush!, - undef_blocks, blocks @recompile_invalidations begin import StateSelection diff --git a/src/structural_transformation/symbolics_tearing.jl b/src/structural_transformation/symbolics_tearing.jl index aefbea0311..f0eb456b46 100644 --- a/src/structural_transformation/symbolics_tearing.jl +++ b/src/structural_transformation/symbolics_tearing.jl @@ -20,6 +20,30 @@ function tearing(sys::AbstractSystem, state = TearingState(sys); mm = nothing, invalidate_cache!(reassemble_alg(state, tearing_result, mm; fully_determined)) end +function safe_isinteger(@nospecialize(x::Number)) + if x isa Int64 + return true + elseif x isa Int32 + return true + elseif x isa BigInt + return typemin(Int) <= x <= typemax(Int) + elseif x isa Float64 + return isinteger(x) && typemin(Int) <= x <= typemax(Int) + elseif x isa Float32 + return isinteger(x) && typemin(Int) <= x <= typemax(Int) + elseif x isa BigFloat + return isinteger(x) && typemin(Int) <= x <= typemax(Int) + elseif x isa Rational{Int64} + return isinteger(x) && typemin(Int) <= x <= typemax(Int) + elseif x isa Rational{Int32} + return isinteger(x) && typemin(Int) <= x <= typemax(Int) + elseif x isa Rational{BigInt} + return isinteger(x) && typemin(Int) <= x <= typemax(Int) + else + return isinteger(x)::Bool && (typemin(Int) <= x)::Bool && (x <= typemax(Int))::Bool + end +end + """ dummy_derivative(sys) @@ -38,10 +62,8 @@ function dummy_derivative(sys, state = TearingState(sys); el = _J[i] Moshi.Match.@match el begin BSImpl.Const(; val) && if val isa Number end => begin - isinteger(val)::Bool || return nothing - val = Int(val) - typemin(Int) <= val <= typemax(Int) || return nothing - J[i] = val + safe_isinteger(val)|| return nothing + J[i] = convert(Int, val)::Int end _ => return nothing end