From 39a334b4b0430ddc32c55e66709af87f1e754062 Mon Sep 17 00:00:00 2001
From: Aayush Sabharwal <aayush.sabharwal@gmail.com>
Date: Fri, 31 Jan 2025 18:11:33 +0530
Subject: [PATCH 1/6] refactor: remove `HomotopyContinuationProblem`
---
.../nonlinear/homotopy_continuation.jl | 90 -------------------
1 file changed, 90 deletions(-)
diff --git a/src/systems/nonlinear/homotopy_continuation.jl b/src/systems/nonlinear/homotopy_continuation.jl
index 09e11199e8..692b0719cd 100644
--- a/src/systems/nonlinear/homotopy_continuation.jl
+++ b/src/systems/nonlinear/homotopy_continuation.jl
@@ -1,93 +1,3 @@
-"""
-$(TYPEDEF)
-
-A type of Nonlinear problem which specializes on polynomial systems and uses
-HomotopyContinuation.jl to solve the system. Requires importing HomotopyContinuation.jl to
-create and solve.
-"""
-struct HomotopyContinuationProblem{uType, H, D, O, SS, U} <:
- SciMLBase.AbstractNonlinearProblem{uType, true}
- """
- The initial values of states in the system. If there are multiple real roots of
- the system, the one closest to this point is returned.
- """
- u0::uType
- """
- A subtype of `HomotopyContinuation.AbstractSystem` to solve. Also contains the
- parameter object.
- """
- homotopy_continuation_system::H
- """
- A function with signature `(u, p) -> resid`. In case of rational functions, this
- is used to rule out roots of the system which would cause the denominator to be
- zero.
- """
- denominator::D
- """
- The `NonlinearSystem` used to create this problem. Used for symbolic indexing.
- """
- sys::NonlinearSystem
- """
- A function which generates and returns observed expressions for the given system.
- """
- obsfn::O
- """
- The HomotopyContinuation.jl solver and start system, obtained through
- `HomotopyContinuation.solver_startsystems`.
- """
- solver_and_starts::SS
- """
- A function which takes a solution of the transformed system, and returns a vector
- of solutions for the original system. This is utilized when converting systems
- to polynomials.
- """
- unpack_solution::U
-end
-
-function HomotopyContinuationProblem(::AbstractSystem, _u0, _p; kwargs...)
- error("HomotopyContinuation.jl is required to create and solve `HomotopyContinuationProblem`s. Please run `Pkg.add(\"HomotopyContinuation\")` to continue.")
-end
-
-"""
- $(TYPEDSIGNATURES)
-
-Utility function for `safe_HomotopyContinuationProblem`, implemented in the extension.
-"""
-function _safe_HomotopyContinuationProblem end
-
-"""
- $(TYPEDSIGNATURES)
-
-Return a `HomotopyContinuationProblem` if the extension is loaded and the system is
-polynomial. If the extension is not loaded, return `nothing`. If the system is not
-polynomial, return the appropriate `NotPolynomialError`.
-"""
-function safe_HomotopyContinuationProblem(sys::NonlinearSystem, args...; kwargs...)
- if Base.get_extension(ModelingToolkit, :MTKHomotopyContinuationExt) === nothing
- return nothing
- end
- return _safe_HomotopyContinuationProblem(sys, args...; kwargs...)
-end
-
-SymbolicIndexingInterface.symbolic_container(p::HomotopyContinuationProblem) = p.sys
-SymbolicIndexingInterface.state_values(p::HomotopyContinuationProblem) = p.u0
-function SymbolicIndexingInterface.set_state!(p::HomotopyContinuationProblem, args...)
- set_state!(p.u0, args...)
-end
-function SymbolicIndexingInterface.parameter_values(p::HomotopyContinuationProblem)
- parameter_values(p.homotopy_continuation_system)
-end
-function SymbolicIndexingInterface.set_parameter!(p::HomotopyContinuationProblem, args...)
- set_parameter!(parameter_values(p), args...)
-end
-function SymbolicIndexingInterface.observed(p::HomotopyContinuationProblem, sym)
- if p.obsfn !== nothing
- return p.obsfn(sym)
- else
- return SymbolicIndexingInterface.observed(p.sys, sym)
- end
-end
-
function contains_variable(x, wrt)
any(y -> occursin(y, x), wrt)
end
From 945999e86fa621ce78a4b29dc66f17871efd357e Mon Sep 17 00:00:00 2001
From: Aayush Sabharwal <aayush.sabharwal@gmail.com>
Date: Fri, 31 Jan 2025 18:14:04 +0530
Subject: [PATCH 2/6] refactor: rearrange code
---
.../nonlinear/homotopy_continuation.jl | 46 +++++++++----------
1 file changed, 23 insertions(+), 23 deletions(-)
diff --git a/src/systems/nonlinear/homotopy_continuation.jl b/src/systems/nonlinear/homotopy_continuation.jl
index 692b0719cd..714f6bec7b 100644
--- a/src/systems/nonlinear/homotopy_continuation.jl
+++ b/src/systems/nonlinear/homotopy_continuation.jl
@@ -45,29 +45,6 @@ function display_reason(reason::NonPolynomialReason.T, sym)
end
end
-"""
- $(TYPEDEF)
-
-Information about an expression about its polynomial nature.
-"""
-mutable struct PolynomialData
- """
- A list of all non-polynomial terms in the expression.
- """
- non_polynomial_terms::Vector{BasicSymbolic}
- """
- Corresponding to `non_polynomial_terms`, a list of reasons why they are
- not polynomial.
- """
- reasons::Vector{NonPolynomialReason.T}
- """
- Whether the polynomial contains parametric exponents of unknowns.
- """
- has_parametric_exponent::Bool
-end
-
-PolynomialData() = PolynomialData(BasicSymbolic[], NonPolynomialReason.T[], false)
-
abstract type PolynomialTransformationError <: Exception end
struct MultivarTerm <: PolynomialTransformationError
@@ -138,6 +115,29 @@ function Base.showerror(io::IO, err::NotPolynomialError)
end
end
+"""
+ $(TYPEDEF)
+
+Information about an expression about its polynomial nature.
+"""
+mutable struct PolynomialData
+ """
+ A list of all non-polynomial terms in the expression.
+ """
+ non_polynomial_terms::Vector{BasicSymbolic}
+ """
+ Corresponding to `non_polynomial_terms`, a list of reasons why they are
+ not polynomial.
+ """
+ reasons::Vector{NonPolynomialReason.T}
+ """
+ Whether the polynomial contains parametric exponents of unknowns.
+ """
+ has_parametric_exponent::Bool
+end
+
+PolynomialData() = PolynomialData(BasicSymbolic[], NonPolynomialReason.T[], false)
+
function is_polynomial!(data, y, wrt)
process_polynomial!(data, y, wrt)
isempty(data.reasons)
From d39727e535489c2a1baa8de67edb5bbda421057e Mon Sep 17 00:00:00 2001
From: Aayush Sabharwal <aayush.sabharwal@gmail.com>
Date: Fri, 31 Jan 2025 18:42:11 +0530
Subject: [PATCH 3/6] fixup! refactor: rearrange code
---
src/systems/nonlinear/homotopy_continuation.jl | 2 +-
1 file changed, 1 insertion(+), 1 deletion(-)
diff --git a/src/systems/nonlinear/homotopy_continuation.jl b/src/systems/nonlinear/homotopy_continuation.jl
index 714f6bec7b..757dfeb1bb 100644
--- a/src/systems/nonlinear/homotopy_continuation.jl
+++ b/src/systems/nonlinear/homotopy_continuation.jl
@@ -96,7 +96,7 @@ end
struct NotPolynomialError <: Exception
transformation_err::Union{PolynomialTransformationError, Nothing}
eq::Vector{Equation}
- data::Vector{PolynomialData}
+ data::Vector{Any} # PolynomialData but then I have to put that struct above this
end
function Base.showerror(io::IO, err::NotPolynomialError)
From cf2d5c658912378d5d48e136100101db35ef1c79 Mon Sep 17 00:00:00 2001
From: Aayush Sabharwal <aayush.sabharwal@gmail.com>
Date: Sat, 1 Feb 2025 00:55:39 +0530
Subject: [PATCH 4/6] feat: record all variables used in polynomial terms in
`PolynomialData`
---
src/systems/nonlinear/homotopy_continuation.jl | 14 +++++++++++---
1 file changed, 11 insertions(+), 3 deletions(-)
diff --git a/src/systems/nonlinear/homotopy_continuation.jl b/src/systems/nonlinear/homotopy_continuation.jl
index 757dfeb1bb..2ef1120d2f 100644
--- a/src/systems/nonlinear/homotopy_continuation.jl
+++ b/src/systems/nonlinear/homotopy_continuation.jl
@@ -121,6 +121,11 @@ end
Information about an expression about its polynomial nature.
"""
mutable struct PolynomialData
+ """
+ A list of the variables in the expression being searched that occur outside of
+ any non-polynomial terms.
+ """
+ solo_terms::Set{BasicSymbolic}
"""
A list of all non-polynomial terms in the expression.
"""
@@ -136,7 +141,7 @@ mutable struct PolynomialData
has_parametric_exponent::Bool
end
-PolynomialData() = PolynomialData(BasicSymbolic[], NonPolynomialReason.T[], false)
+PolynomialData() = PolynomialData(Set{BasicSymbolic}(), BasicSymbolic[], NonPolynomialReason.T[], false)
function is_polynomial!(data, y, wrt)
process_polynomial!(data, y, wrt)
@@ -152,9 +157,12 @@ writing said information to `data`.
function process_polynomial!(data::PolynomialData, x, wrt)
x = unwrap(x)
symbolic_type(x) == NotSymbolic() && return true
- iscall(x) || return true
contains_variable(x, wrt) || return true
- any(isequal(x), wrt) && return true
+ if any(isequal(x), wrt)
+ push!(data.solo_terms, x)
+ return true
+ end
+ iscall(x) || return true
if operation(x) in (*, +, -, /)
# `map` because `all` will early exit, but we want to search
From b790dfbb6a7be172f515f95acd07076dbc4a0d4e Mon Sep 17 00:00:00 2001
From: Aayush Sabharwal <aayush.sabharwal@gmail.com>
Date: Sat, 1 Feb 2025 12:38:08 +0530
Subject: [PATCH 5/6] feat: add significantly improved polynomial
transformation
---
.../nonlinear/homotopy_continuation.jl | 277 +++++++++++-------
1 file changed, 164 insertions(+), 113 deletions(-)
diff --git a/src/systems/nonlinear/homotopy_continuation.jl b/src/systems/nonlinear/homotopy_continuation.jl
index 2ef1120d2f..4b23a742f6 100644
--- a/src/systems/nonlinear/homotopy_continuation.jl
+++ b/src/systems/nonlinear/homotopy_continuation.jl
@@ -47,50 +47,25 @@ end
abstract type PolynomialTransformationError <: Exception end
-struct MultivarTerm <: PolynomialTransformationError
- term::Any
- vars::Any
+struct UnmatchedUnknowns <: PolynomialTransformationError
+ unmatched::Vector{BasicSymbolic}
end
-function Base.showerror(io::IO, err::MultivarTerm)
+function Base.showerror(io::IO, err::UnmatchedUnknowns)
println(io,
- "Cannot convert system to polynomial: Found term $(err.term) which is a function of multiple unknowns $(err.vars).")
+ "Cannot convert system to polynomial: could not find terms to solve for unknowns $(err.unmatched).")
end
-struct MultipleTermsOfSameVar <: PolynomialTransformationError
- terms::Any
- var::Any
+struct UnmatchedTerms <: PolynomialTransformationError
+ unmatched::Vector{BasicSymbolic}
end
-function Base.showerror(io::IO, err::MultipleTermsOfSameVar)
- println(io,
- "Cannot convert system to polynomial: Found multiple non-polynomial terms $(err.terms) involving the same unknown $(err.var).")
-end
-
-struct SymbolicSolveFailure <: PolynomialTransformationError
- term::Any
- var::Any
-end
-
-function Base.showerror(io::IO, err::SymbolicSolveFailure)
- println(io,
- "Cannot convert system to polynomial: Unable to symbolically solve $(err.term) for $(err.var).")
-end
-
-struct NemoNotLoaded <: PolynomialTransformationError end
-
-function Base.showerror(io::IO, err::NemoNotLoaded)
- println(io,
- "ModelingToolkit may be able to solve this system as a polynomial system if `Nemo` is loaded. Run `import Nemo` and try again.")
-end
-
-struct VariablesAsPolyAndNonPoly <: PolynomialTransformationError
- vars::Any
+function Base.showerror(io::IO, err::UnmatchedTerms)
+ println(io, "Cannot convert system to polynomial: too many non-polynomial terms in system. Unmatched terms are $(err.unmatched).")
end
-function Base.showerror(io::IO, err::VariablesAsPolyAndNonPoly)
- println(io,
- "Cannot convert convert system to polynomial: Variables $(err.vars) occur in both polynomial and non-polynomial terms in the system.")
+function no_nemo_warning()
+ @warn "ModelingToolkit may be able to symbolically solve some non-polynomial terms in this system for all roots if `Nemo` is loaded. Run `import Nemo` and try again to enable this functionality and possibly obtain additional roots."
end
struct NotPolynomialError <: Exception
@@ -216,6 +191,24 @@ struct PolynomialTransformationData
inv_term::Vector{BasicSymbolic}
end
+"""
+ $(TYPEDEF)
+
+Information for how to solve for unknowns involved in non-symbolically-solvable
+non-polynomial terms to turn the system into a polynomial. Used in
+`PolynomialTransformation`.
+"""
+struct NonlinearSolveTransformation
+ """
+ The system which solves for the unknowns of the parent system.
+ """
+ sys::NonlinearSystem
+ """
+ The input variables to this system representing solutions of non-polynomial terms.
+ """
+ inputvars::Vector{BasicSymbolic}
+end
+
"""
$(TYPEDEF)
@@ -224,20 +217,19 @@ system.
"""
struct PolynomialTransformation
"""
- Substitutions mapping non-polynomial terms to temporary unknowns. The system
- is a polynomial in the new unknowns. Currently, each non-polynomial term is a
- function of a single unknown of the original system.
+ The stages in which to recover the solution in terms of the original unknowns, in
+ order.
"""
- substitution_rules::Dict{BasicSymbolic, BasicSymbolic}
+ solve_stages::Vector{Any}
"""
- A vector of expressions involving unknowns of the transformed system, mapping
- back to solutions of the original system.
+ The (previous) stages each stage depends on.
"""
- all_solutions::Vector{Vector{BasicSymbolic}}
+ stage_dependencies::Vector{Vector{Int}}
"""
- The new unknowns of the transformed system.
+ Mapping from terms to new unknowns they are replaced by. The system is a
+ polynomial in the new unknowns.
"""
- new_dvs::Vector{BasicSymbolic}
+ substitution_rules::Dict{BasicSymbolic, BasicSymbolic}
"""
The polynomial data for each equation.
"""
@@ -270,76 +262,135 @@ function PolynomialTransformation(sys::NonlinearSystem)
d -> d.non_polynomial_terms, vcat, polydata; init = BasicSymbolic[])
unique!(all_non_poly_terms)
- # each variable can only be replaced by one non-polynomial expression involving
- # that variable. Keep track of this mapping.
- var_to_nonpoly = Dict{BasicSymbolic, PolynomialTransformationData}()
+ all_solo_vars = mapreduce(d -> d.solo_terms, union, polydata; init = Set{BasicSymbolic}())
+
+ # Graph matches variables to candidates for unknowns of the polynomial system that
+ # they occur in. These unknowns can be solo variables that appear outside of
+ # non-polynomial terms in the system, or non-polynomials.
+ graph = BipartiteGraph(length(dvs), 0)
+ # all solo variables are candidates for unknowns
+ graph_srcs = dvs
+ graph_dsts = BasicSymbolic[]
+ for (i, var) in enumerate(dvs)
+ var in all_solo_vars || continue
+ push!(graph_dsts, var)
+ vert = add_vertex!(graph, DST)
+ add_edge!(graph, i, vert)
+ end
- is_poly = true
- transformation_err = nothing
+ # buffer to prevent reallocations
+ dvs_in_term = Set()
+ # for efficient queries
+ dvs_to_src = Dict(graph_srcs .=> eachindex(graph_srcs))
+ # build out graph with other non-polynomial terms
for t in all_non_poly_terms
- # if the term involves multiple unknowns, we can't invert it
- dvs_in_term = map(x -> occursin(x, t), dvs)
- if count(dvs_in_term) > 1
- transformation_err = MultivarTerm(t, dvs[dvs_in_term])
- is_poly = false
- break
- end
- # we already have a substitution solving for `var`
- var = dvs[findfirst(dvs_in_term)]
- if haskey(var_to_nonpoly, var) && !isequal(var_to_nonpoly[var].term, t)
- transformation_err = MultipleTermsOfSameVar([t, var_to_nonpoly[var].term], var)
- is_poly = false
- break
- end
- # we want to solve `term - new_var` for `var`
- new_var = gensym(Symbol(var))
- new_var = unwrap(only(@variables $new_var))
- invterm = Symbolics.ia_solve(
- t - new_var, var; complex_roots = false, periodic_roots = false, warns = false)
- # if we can't invert it, quit
- if invterm === nothing || isempty(invterm)
- transformation_err = SymbolicSolveFailure(t, var)
- is_poly = false
- break
- end
- # `ia_solve` returns lazy terms i.e. `asin(1.0)` instead of `pi/2`
- # this just evaluates the constant expressions
- invterm = Symbolics.substitute.(invterm, (Dict(),))
- # RootsOf implies Symbolics couldn't solve the inner polynomial because
- # `Nemo` wasn't loaded.
- if any(x -> iscall(x) && operation(x) == Symbolics.RootsOf, invterm)
- transformation_err = NemoNotLoaded()
- is_poly = false
- break
+ empty!(dvs_in_term)
+ vars!(dvs_in_term, t)
+ intersect!(dvs_in_term, dvs)
+ push!(graph_dsts, t)
+ vert = add_vertex!(graph, DST)
+ for var in dvs_in_term
+ add_edge!(graph, dvs_to_src[var], vert)
end
- var_to_nonpoly[var] = PolynomialTransformationData(new_var, t, invterm)
end
+ # Match variables to the candidate unknown we'll use to solve for them.
+ # This is a poor man's version of `structural_simplify`, but if we create
+ # and simplify a `NonlinearSystem` it makes doing symbolic solving more
+ # annoying.
+ matching = BipartiteGraphs.complete(maximal_matching(graph))
+ inv_matching = invview(matching)
+ # matching is from destination to source vertices
+ unassigned_dsts = filter(i -> matching[i] == unassigned, š¯‘ vertices(graph))
+ unassigned_srcs = filter(i -> inv_matching[i] == unassigned, š¯‘‘vertices(graph))
+
# return the error instead of throwing it, so the user can choose what to do
# without having to catch the exception
- if !is_poly
- return NotPolynomialError(transformation_err, eqs, polydata)
+ if !isempty(unassigned_srcs)
+ return NotPolynomialError(UnmatchedUnknowns(graph_srcs[unassigned_srcs]), eqs, polydata)
+ end
+ if !isempty(unassigned_dsts)
+ return NotPolynomialError(UnmatchedTerms(graph_dsts[unassigned_dsts]), eqs, polydata)
+ end
+
+ # At this point, the matching is perfect. Find the SCCs so we know
+ # which terms to solve for which variables.
+ digraph = DiCMOBiGraph{false}(graph, matching)
+ var_sccs = Graphs.strongly_connected_components(digraph)
+ foreach(sort!, var_sccs)
+ # construct a condensation graph of the SCCs so we can topologically sort them
+ scc_graph = MatchedCondensationGraph(digraph, var_sccs)
+ toporder = topological_sort(scc_graph)
+ var_sccs = var_sccs[toporder]
+ # get the corresponding terms
+ term_sccs = map(var_sccs) do scc
+ map(scc) do src
+ inv_matching[src]
+ end
end
+ # keep track of which previous SCCs each SCC depends on
+ dependencies = Vector{Int}[]
+ # the method to solve each stage
+ solve_stages = []
+ # mapping from terms to the new unknowns they are replaced by
subrules = Dict{BasicSymbolic, BasicSymbolic}()
- # corresponding to each unknown in `dvs`, the list of its possible solutions
- # in terms of the new unknown.
- combinations = Vector{BasicSymbolic}[]
- new_dvs = BasicSymbolic[]
- for x in dvs
- if haskey(var_to_nonpoly, x)
- _data = var_to_nonpoly[x]
- # map term to new unknown
- subrules[_data.term] = _data.new_var
- push!(combinations, _data.inv_term)
- push!(new_dvs, _data.new_var)
- else
- push!(combinations, BasicSymbolic[x])
- push!(new_dvs, x)
+ # if we've already emitted the no nemo warning
+ warned_no_nemo = false
+ for (i, (vscc, tscc)) in enumerate(zip(var_sccs, term_sccs))
+ # dependencies are simply outneighbors
+ push!(dependencies, collect(Graphs.outneighbors(scc_graph, i)))
+
+ # whether the SCC is solvable with a single variable
+ single_scc_solvable = length(vscc) == 1
+ # for single-variable SCCs, we use `ia_solve`
+ if single_scc_solvable
+ varidx = vscc[]
+ termidx = tscc[]
+ var = graph_srcs[varidx]
+ t = graph_dsts[termidx]
+ # Create a new variable and representing the non-polynomial term...
+ new_var = unwrap(similar_variable(var, Symbol(var)))
+ # ...and solve for `var` in terms of this new variable.
+ invterm = Symbolics.ia_solve(t - new_var, var; complex_roots = false, periodic_roots = false, warns = false)
+ # `ia_solve` returns lazy terms i.e. `asin(1.0)` instead of `pi/2`
+ # this just evaluates the constant expressions
+ invterm = Symbolics.substitute.(invterm, (Dict(),))
+ # if `ia_solve` returns `nothing`, the broadcast above turns it into `(nothing,)`
+ if invterm === (nothing,) || isempty(invterm)
+ # if we can't invert it, quit
+ single_scc_solvable = false
+ elseif any(x -> iscall(x) && operation(x) == Symbolics.RootsOf, invterm)
+ # RootsOf implies Symbolics couldn't solve the inner polynomial because
+ # `Nemo` wasn't loaded.
+ warned_no_nemo || no_nemo_warning()
+ warned_no_nemo = true
+ single_scc_solvable = false
+ else
+ subrules[t] = new_var
+ push!(solve_stages, PolynomialTransformationData(new_var, t, invterm))
+ end
+ end
+
+ # the SCC was solved with a single variable
+ single_scc_solvable && continue
+
+ # Solve using a `NonlinearSolve`.
+ vars = graph_srcs[vscc]
+ ts = graph_dsts[tscc]
+ # the new variables are inputs to the system, so they're parameters
+ new_vars = map(vars) do var
+ toparam.(unwrap.(similar_variable(var, Symbol(var))))
+ end
+ eqs = collect(0 .~ (ts .- new_vars))
+ scc_sys = complete(NonlinearSystem(eqs; name = Symbol(:scc_, i)))
+ push!(solve_stages, NonlinearSolveTransformation(scc_sys, new_vars))
+ for (new_var, t) in zip(new_vars, ts)
+ subrules[t] = new_var
end
end
- all_solutions = vec(collect.(collect(Iterators.product(combinations...))))
- return PolynomialTransformation(subrules, all_solutions, new_dvs, polydata)
+
+ return PolynomialTransformation(solve_stages, dependencies, subrules, polydata)
end
"""
@@ -352,6 +403,15 @@ in the equations, to rule out invalid roots.
struct PolynomialTransformationResult
sys::NonlinearSystem
denominators::Vector{BasicSymbolic}
+ """
+ The stages in which to recover the solution in terms of the original unknowns, in
+ order.
+ """
+ solve_stages::Vector{Any}
+ """
+ The (previous) stages each stage depends on.
+ """
+ stage_dependencies::Vector{Vector{Int}}
end
"""
@@ -367,22 +427,13 @@ function transform_system(sys::NonlinearSystem, transformation::PolynomialTransf
dvs = unknowns(sys)
eqs = full_equations(sys)
polydata = transformation.polydata
- new_dvs = transformation.new_dvs
- all_solutions = transformation.all_solutions
+ new_dvs = collect(values(subrules))
eqs2 = Equation[]
denoms = BasicSymbolic[]
for eq in eqs
t = eq.rhs - eq.lhs
t = Symbolics.fixpoint_sub(t, subrules; maxiters = length(dvs))
- # the substituted variable occurs outside the substituted term
- poly_and_nonpoly = map(dvs) do x
- all(!isequal(x), new_dvs) && occursin(x, t)
- end
- if any(poly_and_nonpoly)
- return NotPolynomialError(
- VariablesAsPolyAndNonPoly(dvs[poly_and_nonpoly]), eqs, polydata)
- end
num, den = handle_rational_polynomials(t, new_dvs; fraction_cancel_fn)
# make factors different elements, otherwise the nonzero factors artificially
# inflate the error of the zero factor.
@@ -403,7 +454,7 @@ function transform_system(sys::NonlinearSystem, transformation::PolynomialTransf
# remove observed equations to avoid adding them in codegen
@set! sys2.observed = Equation[]
@set! sys2.substitutions = nothing
- return PolynomialTransformationResult(sys2, denoms)
+ return PolynomialTransformationResult(sys2, denoms, transformation.solve_stages, transformation.stage_dependencies)
end
"""
From 108b87ab0a53cbf1394907e8abfed9fc6783476b Mon Sep 17 00:00:00 2001
From: Aayush Sabharwal <aayush.sabharwal@gmail.com>
Date: Sat, 1 Feb 2025 12:40:37 +0530
Subject: [PATCH 6/6] refactor: remove `MTKHomotopyContinuationExt`
---
Project.toml | 3 -
ext/MTKHomotopyContinuationExt.jl | 225 -----------------------
src/systems/nonlinear/nonlinearsystem.jl | 9 +-
3 files changed, 1 insertion(+), 236 deletions(-)
delete mode 100644 ext/MTKHomotopyContinuationExt.jl
diff --git a/Project.toml b/Project.toml
index e4b2238170..f086d8c088 100644
--- a/Project.toml
+++ b/Project.toml
@@ -65,7 +65,6 @@ BifurcationKit = "0f109fa4-8a5d-4b75-95aa-f515264e7665"
ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
DeepDiffs = "ab62b9b5-e342-54a8-a765-a90f495de1a6"
FMI = "14a09403-18e3-468f-ad8a-74f8dda2d9ac"
-HomotopyContinuation = "f213a82b-91d6-5c5d-acf7-10f1c761b327"
InfiniteOpt = "20393b10-9daf-11e9-18c9-8db751c92c57"
LabelledArrays = "2ee39098-c373-598a-b85f-a56591580800"
@@ -74,7 +73,6 @@ MTKBifurcationKitExt = "BifurcationKit"
MTKChainRulesCoreExt = "ChainRulesCore"
MTKDeepDiffsExt = "DeepDiffs"
MTKFMIExt = "FMI"
-MTKHomotopyContinuationExt = "HomotopyContinuation"
MTKInfiniteOptExt = "InfiniteOpt"
MTKLabelledArraysExt = "LabelledArrays"
@@ -110,7 +108,6 @@ FunctionWrappers = "1.1"
FunctionWrappersWrappers = "0.1"
FMI = "0.14"
Graphs = "1.5.2"
-HomotopyContinuation = "2.11"
InfiniteOpt = "0.5"
InteractiveUtils = "1"
JuliaFormatter = "1.0.47"
diff --git a/ext/MTKHomotopyContinuationExt.jl b/ext/MTKHomotopyContinuationExt.jl
deleted file mode 100644
index 8f17c05b18..0000000000
--- a/ext/MTKHomotopyContinuationExt.jl
+++ /dev/null
@@ -1,225 +0,0 @@
-module MTKHomotopyContinuationExt
-
-using ModelingToolkit
-using ModelingToolkit.SciMLBase
-using ModelingToolkit.Symbolics: unwrap, symtype, BasicSymbolic, simplify_fractions
-using ModelingToolkit.SymbolicIndexingInterface
-using ModelingToolkit.DocStringExtensions
-using HomotopyContinuation
-using ModelingToolkit: iscomplete, parameters, has_index_cache, get_index_cache, get_u0,
- get_u0_p, check_eqs_u0, CommonSolve
-
-const MTK = ModelingToolkit
-
-"""
-$(TYPEDSIGNATURES)
-
-Convert `expr` from a symbolics expression to one that uses `HomotopyContinuation.ModelKit`.
-"""
-function symbolics_to_hc(expr)
- if iscall(expr)
- if operation(expr) == getindex
- args = arguments(expr)
- return ModelKit.Variable(getname(args[1]), args[2:end]...)
- else
- return operation(expr)(symbolics_to_hc.(arguments(expr))...)
- end
- elseif symbolic_type(expr) == NotSymbolic()
- return expr
- else
- return ModelKit.Variable(getname(expr))
- end
-end
-
-"""
-$(TYPEDEF)
-
-A subtype of `HomotopyContinuation.AbstractSystem` used to solve `HomotopyContinuationProblem`s.
-"""
-struct MTKHomotopySystem{F, P, J, V} <: HomotopyContinuation.AbstractSystem
- """
- The generated function for the residual of the polynomial system. In-place.
- """
- f::F
- """
- The parameter object.
- """
- p::P
- """
- The generated function for the jacobian of the polynomial system. In-place.
- """
- jac::J
- """
- The `HomotopyContinuation.ModelKit.Variable` representation of the unknowns of
- the system.
- """
- vars::V
- """
- The number of polynomials in the system. Must also be equal to `length(vars)`.
- """
- nexprs::Int
-end
-
-Base.size(sys::MTKHomotopySystem) = (sys.nexprs, length(sys.vars))
-ModelKit.variables(sys::MTKHomotopySystem) = sys.vars
-
-function (sys::MTKHomotopySystem)(x, p = nothing)
- sys.f(x, sys.p)
-end
-
-function ModelKit.evaluate!(u, sys::MTKHomotopySystem, x, p = nothing)
- sys.f(u, x, sys.p)
-end
-
-function ModelKit.evaluate_and_jacobian!(u, U, sys::MTKHomotopySystem, x, p = nothing)
- sys.f(u, x, sys.p)
- sys.jac(U, x, sys.p)
-end
-
-SymbolicIndexingInterface.parameter_values(s::MTKHomotopySystem) = s.p
-
-"""
- $(TYPEDSIGNATURES)
-
-Create a `HomotopyContinuationProblem` from a `NonlinearSystem` with polynomial equations.
-The problem will be solved by HomotopyContinuation.jl. The resultant `NonlinearSolution`
-will contain the polynomial root closest to the point specified by `u0map` (if real roots
-exist for the system).
-
-Keyword arguments:
-- `eval_expression`: Whether to `eval` the generated functions or use a `RuntimeGeneratedFunction`.
-- `eval_module`: The module to use for `eval`/`@RuntimeGeneratedFunction`
-- `warn_parametric_exponent`: Whether to warn if the system contains a parametric
- exponent preventing the homotopy from being cached.
-
-All other keyword arguments are forwarded to `HomotopyContinuation.solver_startsystems`.
-"""
-function MTK.HomotopyContinuationProblem(
- sys::NonlinearSystem, u0map, parammap = nothing; kwargs...)
- prob = MTK._safe_HomotopyContinuationProblem(sys, u0map, parammap; kwargs...)
- prob isa MTK.HomotopyContinuationProblem || throw(prob)
- return prob
-end
-
-function MTK._safe_HomotopyContinuationProblem(sys, u0map, parammap = nothing;
- fraction_cancel_fn = SymbolicUtils.simplify_fractions, kwargs...)
- if !iscomplete(sys)
- error("A completed `NonlinearSystem` is required. Call `complete` or `structural_simplify` on the system before creating a `HomotopyContinuationProblem`")
- end
- transformation = MTK.PolynomialTransformation(sys)
- if transformation isa MTK.NotPolynomialError
- return transformation
- end
- result = MTK.transform_system(sys, transformation; fraction_cancel_fn)
- if result isa MTK.NotPolynomialError
- return result
- end
- MTK.HomotopyContinuationProblem(sys, transformation, result, u0map, parammap; kwargs...)
-end
-
-function MTK.HomotopyContinuationProblem(
- sys::MTK.NonlinearSystem, transformation::MTK.PolynomialTransformation,
- result::MTK.PolynomialTransformationResult, u0map,
- parammap = nothing; eval_expression = false,
- eval_module = ModelingToolkit, warn_parametric_exponent = true, kwargs...)
- sys2 = result.sys
- denoms = result.denominators
- polydata = transformation.polydata
- new_dvs = transformation.new_dvs
- all_solutions = transformation.all_solutions
-
- _, u0, p = MTK.process_SciMLProblem(
- MTK.EmptySciMLFunction, sys, u0map, parammap; eval_expression, eval_module)
- nlfn = NonlinearFunction{true}(sys2; jac = true, eval_expression, eval_module)
-
- denominator = MTK.build_explicit_observed_function(sys2, denoms)
- unpack_solution = MTK.build_explicit_observed_function(sys2, all_solutions)
-
- hvars = symbolics_to_hc.(new_dvs)
- mtkhsys = MTKHomotopySystem(nlfn.f, p, nlfn.jac, hvars, length(new_dvs))
-
- obsfn = MTK.ObservedFunctionCache(sys; eval_expression, eval_module)
-
- has_parametric_exponents = any(d -> d.has_parametric_exponent, polydata)
- if has_parametric_exponents
- if warn_parametric_exponent
- @warn """
- The system has parametric exponents, preventing caching of the homotopy. \
- This will cause `solve` to be slower. Pass `warn_parametric_exponent \
- = false` to turn off this warning
- """
- end
- solver_and_starts = nothing
- else
- solver_and_starts = HomotopyContinuation.solver_startsolutions(mtkhsys; kwargs...)
- end
- return MTK.HomotopyContinuationProblem(
- u0, mtkhsys, denominator, sys, obsfn, solver_and_starts, unpack_solution)
-end
-
-"""
-$(TYPEDSIGNATURES)
-
-Solve a `HomotopyContinuationProblem`. Ignores the algorithm passed to it, and always
-uses `HomotopyContinuation.jl`. The original solution as returned by
-`HomotopyContinuation.jl` will be available in the `.original` field of the returned
-`NonlinearSolution`.
-
-All keyword arguments except the ones listed below are forwarded to
-`HomotopyContinuation.solve`. Note that the solver and start solutions are precomputed,
-and only keyword arguments related to the solve process are valid. All keyword
-arguments have their default values in HomotopyContinuation.jl, except `show_progress`
-which defaults to `false`.
-
-Extra keyword arguments:
-- `denominator_abstol`: In case `prob` is solving a rational function, roots which cause
- the denominator to be below `denominator_abstol` will be discarded.
-"""
-function CommonSolve.solve(prob::MTK.HomotopyContinuationProblem,
- alg = nothing; show_progress = false, denominator_abstol = 1e-7, kwargs...)
- if prob.solver_and_starts === nothing
- sol = HomotopyContinuation.solve(
- prob.homotopy_continuation_system; show_progress, kwargs...)
- else
- solver, starts = prob.solver_and_starts
- sol = HomotopyContinuation.solve(solver, starts; show_progress, kwargs...)
- end
- realsols = HomotopyContinuation.results(sol; only_real = true)
- if isempty(realsols)
- u = state_values(prob)
- retcode = SciMLBase.ReturnCode.ConvergenceFailure
- resid = prob.homotopy_continuation_system(u)
- else
- T = eltype(state_values(prob))
- distance = T(Inf)
- u = state_values(prob)
- resid = nothing
- for result in realsols
- if any(<=(denominator_abstol),
- prob.denominator(real.(result.solution), parameter_values(prob)))
- continue
- end
- for truesol in prob.unpack_solution(result.solution, parameter_values(prob))
- dist = norm(truesol - state_values(prob))
- if dist < distance
- distance = dist
- u = T.(real.(truesol))
- resid = T.(real.(prob.homotopy_continuation_system(result.solution)))
- end
- end
- end
- # all roots cause denominator to be zero
- if isinf(distance)
- u = state_values(prob)
- resid = prob.homotopy_continuation_system(u)
- retcode = SciMLBase.ReturnCode.Infeasible
- else
- retcode = SciMLBase.ReturnCode.Success
- end
- end
-
- return SciMLBase.build_solution(
- prob, :HomotopyContinuation, u, resid; retcode, original = sol)
-end
-
-end
diff --git a/src/systems/nonlinear/nonlinearsystem.jl b/src/systems/nonlinear/nonlinearsystem.jl
index 89663f7dbc..3463a9a714 100644
--- a/src/systems/nonlinear/nonlinearsystem.jl
+++ b/src/systems/nonlinear/nonlinearsystem.jl
@@ -527,17 +527,10 @@ end
function DiffEqBase.NonlinearProblem{iip}(sys::NonlinearSystem, u0map,
parammap = DiffEqBase.NullParameters();
- check_length = true, use_homotopy_continuation = false, kwargs...) where {iip}
+ check_length = true, kwargs...) where {iip}
if !iscomplete(sys)
error("A completed `NonlinearSystem` is required. Call `complete` or `structural_simplify` on the system before creating a `NonlinearProblem`")
end
- if use_homotopy_continuation
- prob = safe_HomotopyContinuationProblem(
- sys, u0map, parammap; check_length, kwargs...)
- if prob isa HomotopyContinuationProblem
- return prob
- end
- end
f, u0, p = process_SciMLProblem(NonlinearFunction{iip}, sys, u0map, parammap;
check_length, kwargs...)
pt = something(get_metadata(sys), StandardNonlinearProblem())