Fix HomotopyCtx docs/reexport CI failures on top of #5044

docs/src/basics/Homotopy.md

@@ -59,6 +59,26 @@ keep symbolic jacobians, `tgrad`, and index reduction consistent. At runtime,
 differentiated equations reproduce `actual`'s derivative; along the continuation
 they follow the derivative of the blended expression above.
 
+## Opting Out of the Lowering
+
+Some targets cannot lower to a continuation solver — for example a GPU ensemble
+kernel with a fixed set of solvers. For those, pass `homotopy = false` to
+`mtkcompile`:
+
+```julia
+sys = mtkcompile(sys; homotopy = false)
+```
+
+Every `homotopy(actual, simplified)` node is then replaced by `actual` before
+compilation. The compiled system contains no `homotopy` nodes, so the generated
+code evaluates `actual` directly (the `simplified` expression is never emitted),
+`AbstractNonlinearProblem(sys, op)` and the initialization problem are plain
+nonlinear problems rather than a `HomotopyProblem`, and the initialization and
+event affect systems derived from the compiled system are compiled the same way.
+Only the `simplified` starting heuristic is discarded; the solved system is
+unchanged. `ModelingToolkitBase.homotopy_enabled(sys)` reports the setting a
+compiled system was built with.
+
 ## Building a `HomotopyProblem`
 
 A system whose equations contain `homotopy` nodes is built into a
@@ -152,6 +172,9 @@ defaults.
 
 ```@docs
 ModelingToolkit.homotopy
+ModelingToolkitBase.strip_homotopy
+ModelingToolkitBase.homotopy_enabled
+ModelingToolkitBase.HomotopyCtx
 ModelingToolkitBase.NonPolynomialReason
 ```
 

lib/ModelingToolkitBase/src/ModelingToolkitBase.jl

@@ -511,6 +511,7 @@ const set_scalar_metadata = setmetadata
 @public SymbolicADDisallowed, check_symbolic_ad_allowed
 @public tobrownian, toparam
 @public ProblemTypeCtx
+@public HomotopyCtx, homotopy_enabled, strip_homotopy
 
 for prop in [SYS_PROPS; [:continuous_events, :discrete_events]]
     getter = Symbol(:get_, prop)

lib/ModelingToolkitBase/src/problems/initializationproblem.jl

@@ -200,7 +200,10 @@ function mtkcompile_initialization_system(
         isys::AbstractSystem, sys::AbstractSystem; fully_determined, kwargs...
     )
     try
-        return mtkcompile(isys; fully_determined, split = is_split(sys), kwargs...)
+        return mtkcompile(
+            isys; fully_determined, split = is_split(sys),
+            homotopy = homotopy_enabled(sys), kwargs...
+        )
     catch err
         newerr = with_initialization_context(err, isys, sys; kwargs...)
         newerr === err && rethrow()
@@ -253,7 +256,8 @@ equations an unbalanced initialization system is missing, in the same terms as
 function initialization_system_size(isys::AbstractSystem, sys::AbstractSystem; kwargs...)
     return try
         compiled = mtkcompile(
-            isys; fully_determined = false, split = is_split(sys), kwargs...
+            isys; fully_determined = false, split = is_split(sys),
+            homotopy = homotopy_enabled(sys), kwargs...
         )
         (length(equations(compiled)), length(unknowns(compiled)))
     catch

lib/ModelingToolkitBase/src/systems/callbacks.jl

@@ -219,7 +219,11 @@ function AffectSystem(
     # 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))::System
+    affectsys = (
+        @invokelatest mtkcompile(
+            affectsys; fully_determined = nothing, homotopy = homotopy_enabled(parent_sys)
+        )
+    )::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)

lib/ModelingToolkitBase/src/systems/homotopy_operator.jl

@@ -38,7 +38,10 @@ node independently, and all of them share the single continuation parameter `λ`
 `homotopy(actual, simplified)` stays an opaque symbolic operator through `System`
 construction, `mtkcompile`, and runtime code generation — no continuation
 parameter is injected into the system, and systems that do not use the operator
-go through a byte-identical pipeline. Symbolic differentiation works through the
+go through a byte-identical pipeline. Targets that cannot lower to a continuation
+solver can opt out entirely with `mtkcompile(sys; homotopy = false)`, which
+replaces every node by its `actual` branch before compilation (see
+[`strip_homotopy`](@ref)). Symbolic differentiation works through the
 operator: nodewise derivative rules keep symbolic jacobians, `tgrad`, and index
 reduction consistent — at runtime differentiated equations reproduce `actual`'s
 derivative, and along the continuation they follow the blended expression below.
@@ -166,20 +169,17 @@ function has_any_homotopy(sys::AbstractSystem)
 end
 
 """
-    _rewrite_with_lambda(ir, x, λ)
+    _rewrite_homotopy_nodes(ir, x, replace)
 
 Recursively replace every `homotopy(a, s)` node in the unwrapped expression `x`
-with `(1 - λ)*s + λ*a`, where `λ` is supplied by the caller. At λ=1 the lowered
-expression reduces numerically to `actual` (trivial form); at λ=0 it reduces to
-`simplified`.
+with `replace(a, s)`, walking the dependency graph of `ir` bottom-up so that
+parents of rewritten nodes are rebuilt with the new arguments.
 
 Hand-written TermInterface walk rather than a SymbolicUtils `@rule` — keeps the
 recursion explicit and avoids the rule-rewriter overhead for a single-node rewrite.
 """
-function _rewrite_with_lambda(ir::IRStructure{VartypeT}, x::SymbolicT, λ::SymbolicT)
+function _rewrite_homotopy_nodes(ir::IRStructure{VartypeT}, x::SymbolicT, replace::F) where {F}
     xidx = SU.populate_ir!(ir, x)
-    SU.populate_ir!(ir, λ)
-    oneminus = 1 - λ
     reachability = SU.get_reachability(ir, xidx)
     push!(reachability, xidx)
     old_to_new_idxs = Dict{Int32, Int32}()
@@ -200,7 +200,7 @@ function _rewrite_with_lambda(ir::IRStructure{VartypeT}, x::SymbolicT, λ::Symbo
             args = parent(arguments(sym))
         end
         if op === homotopy
-            new_sym = λ * args[1] + oneminus * args[2]
+            new_sym = replace(args[1], args[2])
             old_to_new_idxs[idx] = SU.populate_ir!(ir, new_sym)
             dirty = true
         elseif dirty
@@ -212,6 +212,99 @@ function _rewrite_with_lambda(ir::IRStructure{VartypeT}, x::SymbolicT, λ::Symbo
     return ir[get(old_to_new_idxs, xidx, xidx)]
 end
 
+"""
+    _rewrite_with_lambda(ir, x, λ)
+
+Recursively replace every `homotopy(a, s)` node in the unwrapped expression `x`
+with `(1 - λ)*s + λ*a`, where `λ` is supplied by the caller. At λ=1 the lowered
+expression reduces numerically to `actual` (trivial form); at λ=0 it reduces to
+`simplified`.
+"""
+function _rewrite_with_lambda(ir::IRStructure{VartypeT}, x::SymbolicT, λ::SymbolicT)
+    SU.populate_ir!(ir, λ)
+    oneminus = 1 - λ
+    return _rewrite_homotopy_nodes(ir, x, (a, s) -> λ * a + oneminus * s)
+end
+
+"""
+    _strip_homotopy(ir, x)
+
+Recursively replace every `homotopy(actual, simplified)` node in the unwrapped
+expression `x` with `actual`, discarding the `simplified` branch. Returns `x`
+itself when it contains no `homotopy` node.
+"""
+function _strip_homotopy(ir::IRStructure{VartypeT}, x::SymbolicT)
+    SU.populate_ir!(ir, x)
+    has_homotopy(ir, x) || return x
+    return _rewrite_homotopy_nodes(ir, x, (a, s) -> a)
+end
+
+"""
+    strip_homotopy(sys::System)
+
+Return a copy of `sys` in which every Modelica `homotopy(actual, simplified)` node
+has been replaced by its `actual` branch, recursively through subsystems. This is
+what `mtkcompile(sys; homotopy = false)` applies before compilation: the returned
+system contains no `homotopy` nodes, so generated code evaluates `actual` directly
+(the `simplified` expression is neither emitted nor evaluated) and problem
+construction never selects a [`SciMLBase.HomotopyProblem`](@ref).
+
+The rewrite covers the equations, observed equations, initialization equations,
+noise equations, costs, constraints and the values of bindings, initial conditions
+and guesses. Symbolic events and jumps are left as-is; a `homotopy` node inside
+them evaluates as `actual` through the numeric fallback.
+"""
+function strip_homotopy(sys::System)
+    ir = get_irstructure(sys)
+    rew(x::SymbolicT) = _strip_homotopy(ir, x)
+    rew(x) = rew(unwrap(x))
+    rew(eq::Equation) = Equation(rew(eq.lhs), rew(eq.rhs))
+    rew(x::Union{Vector, Matrix}) = map(rew, x)
+    rew(x::Nothing) = nothing
+    function rew_map(d::SymmapT)
+        newd = copy(d)
+        for (k, v) in d
+            newv = rew(v)
+            newv === v || (newd.dict[k] = newv)
+        end
+        return newd
+    end
+    @set! sys.eqs = rew(get_eqs(sys))
+    @set! sys.observed = rew(get_observed(sys))
+    @set! sys.initialization_eqs = rew(get_initialization_eqs(sys))
+    @set! sys.noise_eqs = rew(get_noise_eqs(sys))
+    @set! sys.costs = rew(get_costs(sys))
+    @set! sys.constraints = map(get_constraints(sys)) do c
+        c isa Equation ? rew(c) : c
+    end
+    @set! sys.bindings = rew_map(parent(get_bindings(sys)))
+    @set! sys.initial_conditions = rew_map(get_initial_conditions(sys))
+    @set! sys.guesses = rew_map(get_guesses(sys))
+    @set! sys.systems = map(strip_homotopy, get_systems(sys))
+    return sys
+end
+
+"""
+    HomotopyCtx
+
+System metadata key recording the `homotopy` keyword argument that
+[`mtkcompile`](@ref) was called with. `mtkcompile(sys; homotopy = false)` strips
+every `homotopy(actual, simplified)` node down to `actual` and stores `false` under
+this key so that systems derived from the compiled one (the initialization system,
+event affect systems) are compiled the same way. Query it with
+[`homotopy_enabled`](@ref).
+"""
+struct HomotopyCtx end
+
+"""
+    homotopy_enabled(sys::AbstractSystem)
+
+Whether `homotopy(actual, simplified)` nodes in `sys` are lowered to a continuation
+solve. Returns `false` iff `sys` was compiled with `mtkcompile(sys; homotopy = false)`
+(see [`HomotopyCtx`](@ref)).
+"""
+homotopy_enabled(sys::AbstractSystem) = getmetadata(sys, HomotopyCtx, true)::Bool
+
 # The continuation parameter is a fixed sentinel symbol, NOT `gensym`. A
 # `gensym` name embeds a process-global counter, so the same system would lower
 # to a different name (and thus a different generated `Expr`) in the precompile

lib/ModelingToolkitBase/src/systems/systems.jl

@@ -101,6 +101,13 @@ once — calling `mtkcompile` on an already-compiled system throws
 - `split = true`: Whether the compiled system uses the split parameter representation,
   which stores parameters in type-homogeneous buffers indexed by an `IndexCache`. Pass
   `false` to use a flat parameter vector instead.
+- `homotopy = true`: Whether Modelica [`homotopy`](@ref)`(actual, simplified)` operators
+  are kept for lowering to a continuation solve. Pass `false` to replace every such node
+  by its `actual` branch before compilation: the generated code then contains only
+  `actual` (the `simplified` expression is never emitted), problem construction never
+  selects a `SciMLBase.HomotopyProblem`, and the initialization and event affect systems
+  derived from the compiled system are compiled the same way. Use this for targets that
+  cannot lower to a continuation solver. See [`strip_homotopy`](@ref).
 
 Remaining keyword arguments are forwarded to the internal compilation passes.
 
@@ -129,9 +136,12 @@ function mtkcompile(
         sys::System; additional_passes = (),
         inputs = SymbolicT[], outputs = SymbolicT[],
         disturbance_inputs = SymbolicT[],
-        split = true, kwargs...
+        split = true, homotopy = true, kwargs...
     )
     isscheduled(sys) && throw(RepeatedStructuralSimplificationError())
+    if !homotopy
+        sys = strip_homotopy(sys)
+    end
 
     # For backward compatibility with old ModelingToolkit which does not
     # integrate with the reversible transformation API.
@@ -150,6 +160,11 @@ function mtkcompile(
         newsys = pass(newsys)
     end
     @set! newsys.parent = toggle_namespacing(sys, false)
+    # Record the choice so systems derived from `newsys` (initialization system, event
+    # affect systems) are compiled with the same `homotopy` setting.
+    if !homotopy
+        newsys = setmetadata(newsys, HomotopyCtx, false)
+    end
     # Singular systems may end up with parameter-only equations, which shouldn't error on `complete`
     newsys = complete(newsys; split, allow_parameter_eqs = true)
     return newsys

lib/ModelingToolkitBase/test/homotopy_disabled.jl

@@ -0,0 +1,84 @@
+using ModelingToolkitBase
+using ModelingToolkitBase: homotopy, has_homotopy, has_any_homotopy, homotopy_enabled,
+    strip_homotopy, get_nonlinear_problem_type, generate_rhs,
+    t_nounits as t, D_nounits as D
+using SciMLBase
+using Symbolics
+using Test
+
+# `mtkcompile(sys; homotopy = false)` opts out of the continuation lowering: every
+# `homotopy(actual, simplified)` node is replaced by `actual` before compilation, so the
+# generated code only contains `actual`, problem construction never selects a
+# `HomotopyProblem`, and systems derived from the compiled one inherit the setting.
+
+@testset "strip_homotopy rewrites every field, recursively" begin
+    @variables x(t) y(t) z(t)
+    @named inner = System([0 ~ homotopy(y^2 - 2, y - 1.4)], t, [y], [])
+    @named sys = System(
+        [D(x) ~ homotopy(-x, -x + 1), z ~ homotopy(x^2, x)], t;
+        systems = [inner], initialization_eqs = [homotopy(x^2, x) ~ 1],
+        guesses = [x => homotopy(z + 1, z)]
+    )
+    ssys = strip_homotopy(sys)
+    @test !any(eq -> has_homotopy(eq.lhs) || has_homotopy(eq.rhs), equations(ssys))
+    @test !any(eq -> has_homotopy(eq.lhs) || has_homotopy(eq.rhs), initialization_equations(ssys))
+    @test !any(has_homotopy, values(guesses(ssys)))
+    @test isequal(equations(ssys)[1].rhs, -x)
+    @test isequal(equations(ssys)[2].rhs, x^2)
+    @test isequal(initialization_equations(ssys)[1].lhs, x^2)
+    # the original is untouched
+    @test any(eq -> has_homotopy(eq.rhs), equations(sys))
+end
+
+@testset "homotopy = false: generated code uses only the actual branch" begin
+    @variables y
+    @named sys = System([0 ~ homotopy(atan(y - 3), y)])
+    csys = mtkcompile(sys; homotopy = false)
+    @test !homotopy_enabled(csys)
+    @test !has_any_homotopy(csys)
+
+    expr = generate_rhs(csys; expression = Val{true})
+    @test !occursin("homotopy", string(expr))
+    @test occursin("atan", string(expr))
+
+    @test get_nonlinear_problem_type(csys) === NonlinearProblem
+    prob = SciMLBase.AbstractNonlinearProblem(csys, [y => 12.0])
+    @test prob isa NonlinearProblem
+    @test prob.f([12.0], prob.p) ≈ [atan(9.0)]
+    @test_throws ArgumentError HomotopyProblem(csys, [y => 12.0])
+
+    # default: nodes are kept and lowered
+    dsys = mtkcompile(sys)
+    @test homotopy_enabled(dsys)
+    @test has_any_homotopy(dsys)
+    @test occursin("homotopy", string(generate_rhs(dsys; expression = Val{true})))
+    @test SciMLBase.AbstractNonlinearProblem(dsys, [y => 12.0]) isa HomotopyProblem
+end
+
+@testset "homotopy = false: initialization inherits the setting" begin
+    @variables x(t) y(t)
+    @named sys = System([D(x) ~ -x, 0 ~ homotopy(y^2 - x, y - 1)], t; guesses = [y => 1.5])
+    csys = mtkcompile(sys; homotopy = false)
+    prob = ODEProblem(csys, [x => 1.0], (0.0, 1.0))
+    iprob = prob.f.initialization_data.initializeprob
+    @test iprob isa NonlinearProblem
+    @test !occursin("homotopy", string(generate_rhs(iprob.f.sys; expression = Val{true})))
+
+    # initialization equations supplied at problem construction are stripped as well
+    @named sys2 = System([D(x) ~ -x], t)
+    csys2 = mtkcompile(sys2; homotopy = false)
+    prob2 = ODEProblem(
+        csys2, [], (0.0, 1.0);
+        initialization_eqs = [homotopy(atan(x - 3) + x - 3, x - 3) ~ 0], guesses = [x => 1.0]
+    )
+    iprob2 = prob2.f.initialization_data.initializeprob
+    @test iprob2 isa NonlinearProblem
+    @test !has_any_homotopy(iprob2.f.sys)
+
+    # with the default, the same problem routes to a HomotopyProblem
+    prob3 = ODEProblem(
+        mtkcompile(sys2), [], (0.0, 1.0);
+        initialization_eqs = [homotopy(atan(x - 3) + x - 3, x - 3) ~ 0], guesses = [x => 1.0]
+    )
+    @test prob3.f.initialization_data.initializeprob isa HomotopyProblem
+end

lib/ModelingToolkitBase/test/runtests.jl

@@ -98,6 +98,7 @@ end
         @safetestset "Homotopy problem construction & sweep" include("homotopy_problem.jl")
         @safetestset "Homotopy OMC parity" include("homotopy_omc_parity.jl")
         @safetestset "Homotopy initialization routing" include("homotopy_initialization.jl")
+        @safetestset "Homotopy disabled in mtkcompile" include("homotopy_disabled.jl")
         @safetestset "PDE Construction Test" include("pdesystem.jl")
         @safetestset "JumpSystem Test" include("jumpsystem.jl")
         @safetestset "Poissonians Test" include("poissonians.jl")

test/homotopy_initialization_scc.jl

@@ -117,3 +117,23 @@ end
     @test SciMLBase.successful_retcode(sol)
     @test sol[a] ≈ sqrt(2) atol = 1.0e-8
 end
+
+@testset "homotopy = false: SCC init has no HomotopyProblem blocks" begin
+    @variables x(t) a(t) b(t) c(t)
+    @named sys = System(
+        [D(x) ~ -x, 0 ~ homotopy(a^2 - 2, a - 1.414), 0 ~ b^2 - a, 0 ~ c^2 - b],
+        t; guesses = [a => 1.5, b => 1.2, c => 1.1]
+    )
+    sys = mtkcompile(sys; homotopy = false)
+    @test !ModelingToolkit.homotopy_enabled(sys)
+    @test !ModelingToolkit.has_any_homotopy(sys)
+    iprob = ModelingToolkit.InitializationProblem(
+        sys, 0.0, [x => 1.0]; warn_initialize_determined = false
+    )
+    @test iprob isa SciMLBase.SCCNonlinearProblem
+    @test !any(p -> p isa SciMLBase.HomotopyProblem, iprob.probs)
+    sol = solve(iprob)
+    @test SciMLBase.successful_retcode(sol)
+    @test sol[a] ≈ sqrt(2) atol = 1.0e-6
+    @test sol[c] ≈ sqrt(sqrt(sqrt(2))) atol = 1.0e-6
+end