Skip to content

Commit 07aa824

Browse files
docs: add documentation for HomotopyContinuation interface
1 parent 6667457 commit 07aa824

File tree

2 files changed

+83
-9
lines changed

2 files changed

+83
-9
lines changed

ext/MTKHomotopyContinuationExt.jl

Lines changed: 58 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -4,18 +4,22 @@ using ModelingToolkit
44
using ModelingToolkit.SciMLBase
55
using ModelingToolkit.Symbolics: unwrap
66
using ModelingToolkit.SymbolicIndexingInterface
7+
using ModelingToolkit.DocStringExtensions
78
using HomotopyContinuation
89
using ModelingToolkit: iscomplete, parameters, has_index_cache, get_index_cache, get_u0,
910
get_u0_p, check_eqs_u0, CommonSolve
1011

1112
const MTK = ModelingToolkit
1213

1314
function contains_variable(x, wrt)
14-
any(isequal(x), wrt) && return true
15-
iscall(x) || return false
16-
return any(y -> contains_variable(y, wrt), arguments(x))
15+
any(y -> occursin(y, x), wrt)
1716
end
1817

18+
"""
19+
$(TYPEDSIGNATURES)
20+
21+
Check if `x` is polynomial with respect to the variables in `wrt`.
22+
"""
1923
function is_polynomial(x, wrt)
2024
x = unwrap(x)
2125
symbolic_type(x) == NotSymbolic() && return true
@@ -33,6 +37,11 @@ function is_polynomial(x, wrt)
3337
return false
3438
end
3539

40+
"""
41+
$(TYPEDSIGNATURES)
42+
43+
Convert `expr` from a symbolics expression to one that uses `HomotopyContinuation.ModelKit`.
44+
"""
3645
function symbolics_to_hc(expr)
3746
if iscall(expr)
3847
if operation(expr) == getindex
@@ -48,11 +57,32 @@ function symbolics_to_hc(expr)
4857
end
4958
end
5059

60+
"""
61+
$(TYPEDEF)
62+
63+
A subtype of `HomotopyContinuation.AbstractSystem` used to solve `HomotopyContinuationProblem`s.
64+
"""
5165
struct MTKHomotopySystem{F, P, J, V} <: HomotopyContinuation.AbstractSystem
66+
"""
67+
The generated function for the residual of the polynomial system. In-place.
68+
"""
5269
f::F
70+
"""
71+
The parameter object.
72+
"""
5373
p::P
74+
"""
75+
The generated function for the jacobian of the polynomial system. In-place.
76+
"""
5477
jac::J
78+
"""
79+
The `HomotopyContinuation.ModelKit.Variable` representation of the unknowns of
80+
the system.
81+
"""
5582
vars::V
83+
"""
84+
The number of polynomials in the system. Must also be equal to `length(vars)`.
85+
"""
5686
nexprs::Int
5787
end
5888

@@ -74,8 +104,17 @@ end
74104

75105
SymbolicIndexingInterface.parameter_values(s::MTKHomotopySystem) = s.p
76106

107+
"""
108+
$(TYPEDSIGNATURES)
109+
110+
Create a `HomotopyContinuationProblem` from a `NonlinearSystem` with polynomial equations.
111+
The problem will be solved by HomotopyContinuation.jl. The resultant `NonlinearSolution`
112+
will contain the polynomial root closest to the point specified by `u0map` (if real roots
113+
exist for the system).
114+
"""
77115
function MTK.HomotopyContinuationProblem(
78-
sys::NonlinearSystem, u0map, parammap; compile = :all, eval_expression = false, eval_module = ModelingToolkit, kwargs...)
116+
sys::NonlinearSystem, u0map, parammap = nothing; eval_expression = false,
117+
eval_module = ModelingToolkit, kwargs...)
79118
if !iscomplete(sys)
80119
error("A completed `NonlinearSystem` is required. Call `complete` or `structural_simplify` on the system before creating a `HomotopyContinuationProblem`")
81120
end
@@ -100,8 +139,21 @@ function MTK.HomotopyContinuationProblem(
100139
return MTK.HomotopyContinuationProblem(u0, mtkhsys, sys, obsfn)
101140
end
102141

103-
function CommonSolve.solve(prob::MTK.HomotopyContinuationProblem; kwargs...)
104-
sol = HomotopyContinuation.solve(prob.homotopy_continuation_system; kwargs...)
142+
"""
143+
$(TYPEDSIGNATURES)
144+
145+
Solve a `HomotopyContinuationProblem`. Ignores the algorithm passed to it, and always
146+
uses `HomotopyContinuation.jl`. All keyword arguments are forwarded to
147+
`HomotopyContinuation.solve`. The original solution as returned by `HomotopyContinuation.jl`
148+
will be available in the `.original` field of the returned `NonlinearSolution`.
149+
150+
All keyword arguments have their default values in HomotopyContinuation.jl, except
151+
`show_progress` which defaults to `false`.
152+
"""
153+
function CommonSolve.solve(prob::MTK.HomotopyContinuationProblem,
154+
alg = nothing; show_progress = false, kwargs...)
155+
sol = HomotopyContinuation.solve(
156+
prob.homotopy_continuation_system; show_progress, kwargs...)
105157
realsols = HomotopyContinuation.results(sol; only_real = true)
106158
if isempty(realsols)
107159
u = state_values(prob)

src/systems/nonlinear/nonlinearsystem.jl

Lines changed: 25 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -566,15 +566,37 @@ function Base.:(==)(sys1::NonlinearSystem, sys2::NonlinearSystem)
566566
all(s1 == s2 for (s1, s2) in zip(get_systems(sys1), get_systems(sys2)))
567567
end
568568

569-
struct HomotopyContinuationProblem{uType, H, O} <: SciMLBase.AbstractNonlinearProblem{uType, true}
569+
"""
570+
$(TYPEDEF)
571+
572+
A type of Nonlinear problem which specializes on polynomial systems and uses
573+
HomotopyContinuation.jl to solve the system. Requires importing HomotopyContinuation.jl to
574+
create and solve.
575+
"""
576+
struct HomotopyContinuationProblem{uType, H, O} <:
577+
SciMLBase.AbstractNonlinearProblem{uType, true}
578+
"""
579+
The initial values of states in the system. If there are multiple real roots of
580+
the system, the one closest to this point is returned.
581+
"""
570582
u0::uType
583+
"""
584+
A subtype of `HomotopyContinuation.AbstractSystem` to solve. Also contains the
585+
parameter object.
586+
"""
571587
homotopy_continuation_system::H
588+
"""
589+
The `NonlinearSystem` used to create this problem. Used for symbolic indexing.
590+
"""
572591
sys::NonlinearSystem
592+
"""
593+
A function which generates and returns observed expressions for the given system.
594+
"""
573595
obsfn::O
574596
end
575597

576-
function HomotopyContinuationProblem(args...; kwargs...)
577-
error("Requires HomotopyContinuationExt")
598+
function HomotopyContinuationProblem(::AbstractSystem, _u0, _p; kwargs...)
599+
error("HomotopyContinuation.jl is required to create and solve `HomotopyContinuationProblem`s. Please run `Pkg.add(\"HomotopyContinuation\")` to continue.")
578600
end
579601

580602
SymbolicIndexingInterface.symbolic_container(p::HomotopyContinuationProblem) = p.sys

0 commit comments

Comments
 (0)