Make simplification deterministic w.r.t. hash/objectid ordering

lib/ModelingToolkitBase/src/ModelingToolkitBase.jl

@@ -229,6 +229,7 @@ include("systems/pde/pdesystem.jl")
 include("systems/unit_check.jl")
 include("systems/dependency_graphs.jl")
 include("discretedomain.jl")
+include("canonical_ordering.jl")
 include("systems/systems.jl")
 
 include("debugging.jl")

lib/ModelingToolkitBase/src/canonical_ordering.jl

@@ -0,0 +1,122 @@
+"""
+    $TYPEDSIGNATURES
+
+If `x` is a small integer (fits in `Int8`) return `Int(x)`. Otherwise, return
+`typemin(Int)`. Used to build the integer index portion of [`canonical_sort_key`](@ref)
+without allocating.
+"""
+function manual_dispatch_is_small_int(@nospecialize(x::Number))::Int
+    if x isa Int
+        return typemin(Int8) <= x <= typemax(Int8) ? x : typemin(Int)
+    elseif x isa BigInt
+        return typemin(Int8) <= x <= typemax(Int8) ? Int(x) : typemin(Int)
+    elseif x isa Int32
+        return typemin(Int8) <= x <= typemax(Int8) ? Int(x) : typemin(Int)
+    elseif x isa Float64
+        return isinteger(x) && typemin(Int8) <= x <= typemax(Int8) ? Int(x) : typemin(Int)
+    elseif x isa Float32
+        return isinteger(x) && typemin(Int8) <= x <= typemax(Int8) ? Int(x) : typemin(Int)
+    elseif x isa BigFloat
+        return isinteger(x) && typemin(Int8) <= x <= typemax(Int8) ? Int(x) : typemin(Int)
+    elseif x isa Rational{Int}
+        return isinteger(x) && typemin(Int8) <= x <= typemax(Int8) ? Int(x) : typemin(Int)
+    elseif x isa Rational{Int32}
+        return isinteger(x) && typemin(Int8) <= x <= typemax(Int8) ? Int(x) : typemin(Int)
+    elseif x isa Rational{BigInt}
+        return isinteger(x) && typemin(Int8) <= x <= typemax(Int8) ? Int(x) : typemin(Int)
+    else
+        return if isinteger(x)::Bool && (typemin(Int8) <= x <= typemax(Int8))::Bool
+            Int(x)::Int
+        else
+            typemin(Int)
+        end
+    end
+end
+
+"""
+    $TYPEDSIGNATURES
+
+Total, allocation-light "canonical name" for any expression that may appear as a
+variable/parameter: the name for named variables (`Sym` / call-variable / `getindex`),
+the operation's name for operator/function applications, and a fixed sentinel symbol for
+the remaining structural variants. Key collisions are acceptable — they only mean the
+canonical tie-break falls back to the original order among the colliding entries.
+"""
+function canonical_name(x::SymbolicT)
+    Moshi.Match.@match x begin
+        BSImpl.Sym(; name) => name
+        BSImpl.Term(; f, args) && if f === getindex end => canonical_name(args[1])
+        BSImpl.Term(; f) => canonical_opname(f)
+        BSImpl.AddMul(; variant) => variant === SU.AddMulVariant.ADD ? :+ : :*
+        BSImpl.Div(;) => :/
+        BSImpl.ArrayOp(;) => Symbol("#arrayop")
+        BSImpl.Const(;) => Symbol("#const")
+        _ => Symbol("#expr")
+    end
+end
+
+function canonical_opname(@nospecialize(f))
+    f isa SymbolicT && return canonical_name(f)
+    f isa Function && return nameof(f)::Symbol
+    return nameof(typeof(f))::Symbol
+end
+
+"""
+    $TYPEDSIGNATURES
+
+Structured, hash-independent canonical sort key for a symbolic variable/parameter: a tuple
+`(name, indices, opsig)` where `name` is the [`canonical_name`](@ref) of the underlying
+(array) variable, `indices` are the integer indices when `v` is a scalarized array element
+(empty otherwise) and `opsig` encodes the operator chain wrapping the variable
+(`Differential` → `1`; `Shift` → `2` followed by its step count; any other single-argument
+operator → `3`). Comparing these tuples orders variables deterministically regardless of
+hashing/`objectid`, declaration order or equation order, without stringifying symbolic
+expressions. Use this (rather than `Dict`/`Set` iteration order or `hash`) wherever the
+order of a collection of symbolic quantities reaches observable output.
+"""
+function canonical_sort_key(v::SymbolicT)
+    # `opsig`/`idxs` are built as `Vector{Int}` (not growing tuples) so the key type is
+    # concrete and inferrable after the loop. Vectors compare lexicographically, so the
+    # tuple ordering is unchanged.
+    opsig = Int[]
+    x = v
+    while true
+        stripped = Moshi.Match.@match x begin
+            BSImpl.Term(; f, args) && if f isa Differential end => begin
+                push!(opsig, 1)
+                args[1]
+            end
+            BSImpl.Term(; f, args) && if f isa Shift end => begin
+                push!(opsig, 2, Int(f.steps))
+                args[1]
+            end
+            BSImpl.Term(; f, args) && if f isa SU.Operator && length(args) == 1 end => begin
+                push!(opsig, 3)
+                args[1]
+            end
+            _ => nothing
+        end
+        stripped === nothing && break
+        x = stripped
+    end
+    idxs = Int[]
+    Moshi.Match.@match x begin
+        BSImpl.Term(; f, args) && if f === getindex end => begin
+            for i in Iterators.drop(args, 1)
+                ival = SU.isconst(i) ? manual_dispatch_is_small_int(unwrap_const(i))::Int : 0
+                push!(idxs, ival)
+            end
+            x = args[1]
+        end
+        _ => nothing
+    end
+    return (canonical_name(x), idxs, opsig)
+end
+
+"""
+    $TYPEDSIGNATURES
+
+Return a new vector containing the symbolic quantities `xs` in deterministic
+[`canonical_sort_key`](@ref) order. Stable: ties retain their original relative order.
+"""
+canonical_sort(xs) = sort(collect(xs); by = canonical_sort_key, alg = Base.Sort.DEFAULT_STABLE)

lib/ModelingToolkitBase/src/systems/index_cache.jl

@@ -125,15 +125,25 @@ function IndexCache(sys::AbstractSystem)
     events = Union{SymbolicContinuousCallback, SymbolicDiscreteCallback}[cevs; devs]
     parse_callbacks_for_discretes!(sys, cevs, disc_param_callbacks, constant_buffers, nonnumeric_buffers, 0)
     parse_callbacks_for_discretes!(sys, devs, disc_param_callbacks, constant_buffers, nonnumeric_buffers, length(cevs))
-    clock_partitions = unique(collect(values(disc_param_callbacks)))::Vector{BitSet}
-    disc_symtypes = Set{TypeT}()
-    for x in keys(disc_param_callbacks)
-        push!(disc_symtypes, symtype(x))
+    # The discrete-parameter buffer layout must be deterministic. `Dict`/`Set` iteration
+    # over the discrete parameters (`keys`) or their symtypes (`Set{TypeT}`) is
+    # objectid-based and unstable for freshly-created-per-compile parameter types, so
+    # order everything by the stable `canonical_sort_key` of the parameters: partitions by
+    # their sorted contents, symtypes by first appearance in canonical parameter order, and
+    # the parameters within each symtype likewise.
+    clock_partitions = sort!(unique(collect(values(disc_param_callbacks))); by = collect)::Vector{BitSet}
+    sorted_disc_params = canonical_sort(keys(disc_param_callbacks))
+    disc_symtypes = TypeT[]
+    disc_symtype_idx = Dict{TypeT, Int}()
+    for sym in sorted_disc_params
+        st = symtype(sym)::TypeT
+        if !haskey(disc_symtype_idx, st)
+            push!(disc_symtypes, st)
+            disc_symtype_idx[st] = length(disc_symtypes)
+        end
     end
-    disc_symtypes = collect(disc_symtypes)::Vector{TypeT}
-    disc_symtype_idx = Dict{TypeT, Int}(zip(disc_symtypes, eachindex(disc_symtypes)))
     disc_syms_by_symtype = [SymbolicParam[] for _ in disc_symtypes]
-    for sym in keys(disc_param_callbacks)
+    for sym in sorted_disc_params
         push!(disc_syms_by_symtype[disc_symtype_idx[symtype(sym)]], sym)
     end
     disc_syms_by_symtype_by_partition = [Vector{SymbolicParam}[] for _ in disc_symtypes]
@@ -468,7 +478,16 @@ end
 function get_buffer_sizes_and_idxs(::Type{BufT}, sys::AbstractSystem, buffers::Dict) where {BufT}
     idxs = BufT()
     buffer_sizes = BufferTemplate[]
-    for (i, (T, buf)) in enumerate(buffers)
+    # The buffer layout (which type maps to which buffer, and which symbol to which slot)
+    # must be deterministic. Iterating `buffers::Dict{TypeT,Set}` directly is NOT: the buffer
+    # types can be freshly-created per compile (e.g. `FunctionWrapper`/closure types), so
+    # their `Dict`/`Set` iteration order is `objectid`-based and varies run-to-run within a
+    # single session. Order the buckets (and the symbols within each) by the stable
+    # `canonical_sort_key` of the symbolic parameters instead.
+    sorted_buffers = [(T, canonical_sort(buf)) for (T, buf) in buffers]
+    sort!(sorted_buffers; by = b -> canonical_sort_key(first(b[2])),
+        alg = Base.Sort.DEFAULT_STABLE)
+    for (i, (T, buf)) in enumerate(sorted_buffers)
         for (j, p) in enumerate(buf)
             ttp = default_toterm(p)
             rp = renamespace(sys, p)

src/problems/sccnonlinearproblem.jl

@@ -377,7 +377,10 @@ function build_caches!(sys::System, decomposition::SCCDecomposition)
         push!(decomposition.scc_cachevars, cachevars)
         push!(decomposition.scc_cacheexprs, cacheexprs)
 
-        for (k, v) in state
+        # Iterate `state` in canonical (hash/objectid-independent) order so the cache
+        # buffers, the generated cache-writer code and `cachetypes` are deterministic
+        # rather than ordered by the `Dict` hash of the (symbolic) keys / (type) buckets.
+        for (k, v) in sort!(collect(state); by = p -> MTKBase.canonical_sort_key(unwrap(first(p))))
             k = unwrap(k)
             v = unwrap(v)
             T = symtype(k)
@@ -386,9 +389,12 @@ function build_caches!(sys::System, decomposition::SCCDecomposition)
             buf = get!(() -> SymbolicT[], cacheexprs, T)
             push!(buf, k)
         end
-        all_cacheexprs = reduce(vcat, values(cacheexprs); init = SymbolicT[])
+        sorted_cache_types = sort!(collect(keys(cachevars)); by = string)
+        all_cacheexprs = reduce(
+            vcat, (cacheexprs[T] for T in sorted_cache_types); init = SymbolicT[])
         # update the sizes of cache buffers
-        for (T, buf) in cachevars
+        for T in sorted_cache_types
+            buf = cachevars[T]
             idx = findfirst(isequal(T), decomposition.cachetypes)
             if idx === nothing
                 push!(decomposition.cachetypes, T)

src/systems/if_lifting.jl

@@ -17,11 +17,14 @@ struct CondRewriter
     is `true`). `expression < 0` is evaluated on an up-crossing and `expression <= 0` is
     evaluated on a down-crossing to get the updated value of the condition variable.
     """
-    conditions::Dict{Any, @NamedTuple{dependency, expression}}
+    conditions::OrderedCollections.OrderedDict{Any, @NamedTuple{dependency, expression}}
 end
 
 function CondRewriter(iv)
-    return CondRewriter(iv, Dict())
+    # Insertion-ordered so that the generated condition variables, callbacks and
+    # parameters are emitted in a deterministic (creation) order rather than the
+    # hash-dependent order of a plain `Dict`.
+    return CondRewriter(iv, OrderedCollections.OrderedDict())
 end
 
 """
@@ -43,8 +46,10 @@ function new_cond_sym(cw::CondRewriter, expr, dep)
         cw.conditions[existing_var] = (dependency = (dep | existing_dep), expression = expr)
         return existing_var
     end
-    # generate a new condition variable
-    cvar = gensym("cond")
+    # generate a new condition variable. Use a deterministic, insertion-order-based name
+    # (rather than `gensym`, whose global counter makes the generated variable names — and
+    # thus the simplified system — differ run-to-run).
+    cvar = Symbol(:ifelse_cond_, length(cw.conditions) + 1)
     st = symtype(expr)
     iv = cw.iv
     cv = unwrap(first(@discretes $(cvar)(iv)::st = true)) # TODO: real init

test/group_interfaceii.jl

@@ -9,6 +9,7 @@ include("shared/mtktestset.jl")
     @mtktestset("NonlinearSystem Test", "nonlinearsystem.jl")
     @safetestset "SCCNonlinearProblem Test" include("scc_nonlinear_problem.jl")
     @safetestset "IfLifting Test" include("if_lifting.jl")
+    @safetestset "Simplification determinism" include("nondeterminism_simplification.jl")
     @mtktestset("Analysis Points Test", "analysis_points.jl")
     @mtktestset("Causal Variables Connection Test", "causal_variables_connection.jl")
     @safetestset "Subsystem replacement" include("substitute_component.jl")

test/nondeterminism_simplification.jl

@@ -0,0 +1,81 @@
+using ModelingToolkit
+using ModelingToolkit: t_nounits as t, D_nounits as D, IfLifting, get_index_cache
+using Test
+import Symbolics
+
+# Regression tests for hash/objectid-order non-determinism in simplification.
+#
+# Symbolic hashing (and `objectid` of freshly-created-per-compile types) is not stable
+# across sessions / SymbolicUtils versions / GC of the weak hash-cons cache. Several
+# pipeline stages used to order observable output (parameter buffer layout, generated
+# callbacks/parameters) by iterating `Dict`/`Set` collections keyed by such values, so the
+# compiled structure could differ run-to-run within a single session. These tests compile
+# the same model many times (with `GC.gc(true)` between runs to provoke hash-cons eviction)
+# and assert the compiled structure is identical every time.
+
+const NRUNS = 150
+
+"""Compile `build()` `NRUNS` times with GC stress; return the number of distinct
+`dump`-strings observed (should be 1)."""
+function ndistinct(build, dump)
+    seen = Set{UInt}()
+    for _ in 1:NRUNS
+        push!(seen, hash(dump(build())))
+        GC.gc(true)
+    end
+    return length(seen)
+end
+
+@testset "Deterministic nonnumeric/constant parameter buffer layout (index_cache)" begin
+    # Parameters whose symtypes are freshly created on every build mimic the
+    # FunctionWrapper/closure parameter types that multibody-style models route into the
+    # nonnumeric buffers. The `Dict{Type,Set}` that lays out those buffers must not be
+    # iterated in `objectid` order.
+    function build()
+        @variables x(t) = 1.0
+        @parameters k = 1.0
+        ps = ModelingToolkit.SymbolicParam[k]
+        for i in 1:6
+            # a genuinely fresh, empty struct type on every build
+            nm = gensym(:NN)
+            FT = @eval Main (struct $nm end; $nm)
+            p = ModelingToolkit.toparam(Symbolics.unwrap(Symbolics.variable(Symbol(:nn, i); T = FT)))
+            push!(ps, ModelingToolkit.setdefault(p, FT.instance))
+        end
+        return mtkcompile(System([D(x) ~ -k * x], t, [x], ps; name = :sys))
+    end
+    function dump(ss)
+        ic = get_index_cache(ss)
+        io = IOBuffer()
+        for p in sort(parameters(ss); by = string)
+            println(io, p, " => ", ModelingToolkit.parameter_index(ss, p))
+        end
+        # normalize the fresh gensym type names so only ORDER differences remain
+        replace(String(take!(io)), r"##NN#\d+" => "NN")
+    end
+    @test ndistinct(build, dump) == 1
+end
+
+@testset "Deterministic if-lifting (condition vars, parameters, callbacks)" begin
+    function build()
+        @variables x1(t)=0.1 x2(t)=0.2 x3(t)=0.3 x4(t)=0.4
+        @parameters a=1.0 b=2.0 c=3.0
+        eqs = [
+            D(x1) ~ ifelse(x2 < 0.5, a * x2, -a * x2) + ifelse(x3 > 0.1, 0.1, -0.1)
+            D(x2) ~ ifelse(x1 < 0.0, b, -b) + ifelse(x4 < x3, 1.0, 2.0)
+            D(x3) ~ ifelse(x2 < x1, c * x1, c * x4) - ifelse(x4 > 0.2, x3, -x3)
+            D(x4) ~ ifelse(x3 < 0.3, x1, x2) + ifelse(x1 > x4, 0.5, -0.5)
+        ]
+        return mtkcompile(System(eqs, t; name = :top); additional_passes = [IfLifting])
+    end
+    function dump(ss)
+        io = IOBuffer()
+        for sec in (equations, unknowns, parameters, observed)
+            for e in sec(ss)
+                println(io, string(e))
+            end
+        end
+        String(take!(io))
+    end
+    @test ndistinct(build, dump) == 1
+end