@@ -4,18 +4,22 @@ using ModelingToolkit
44using ModelingToolkit. SciMLBase
55using ModelingToolkit. Symbolics: unwrap
66using ModelingToolkit. SymbolicIndexingInterface
7+ using ModelingToolkit. DocStringExtensions
78using HomotopyContinuation
89using ModelingToolkit: iscomplete, parameters, has_index_cache, get_index_cache, get_u0,
910 get_u0_p, check_eqs_u0, CommonSolve
1011
1112const MTK = ModelingToolkit
1213
1314function 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)
1716end
1817
18+ """
19+ $(TYPEDSIGNATURES)
20+
21+ Check if `x` is polynomial with respect to the variables in `wrt`.
22+ """
1923function 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
3438end
3539
40+ """
41+ $(TYPEDSIGNATURES)
42+
43+ Convert `expr` from a symbolics expression to one that uses `HomotopyContinuation.ModelKit`.
44+ """
3645function symbolics_to_hc (expr)
3746 if iscall (expr)
3847 if operation (expr) == getindex
@@ -48,11 +57,32 @@ function symbolics_to_hc(expr)
4857 end
4958end
5059
60+ """
61+ $(TYPEDEF)
62+
63+ A subtype of `HomotopyContinuation.AbstractSystem` used to solve `HomotopyContinuationProblem`s.
64+ """
5165struct 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
5787end
5888
74104
75105SymbolicIndexingInterface. 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+ """
77115function 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)
101140end
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)
0 commit comments