juliahomotopycontinuation/homotopycontinuation.jl

A Julia package for solving systems of polynomials via homotopy continuation.

HomotopyContinuation.jl is a Julia package for solving systems of polynomial equations by numerical homotopy continuation.

---

See juliahomotopycontinuation.org for installation instructions, full content overview and detailed documentation.

---

Basic usage

HomotopyContinuation.jl aims at having easy-to-understand top-level commands. Here is a simple example:

using HomotopyContinuation
@var x y; # declare the variables x and y
F = System([x^2+2y, y^2-2])
result = solve(F)

Result with 4 solutions
==================================
• 4 non-singular solutions (2 real)
• 0 singular solutions (0 real)
• 4 paths tracked
• random seed: 902575

For more see our user guides.

Citing HomotopyContinuation.jl

If you find HomotopyContinuation.jl useful in your work, we kindly request that you cite the following paper:

@inproceedings{HomotopyContinuation.jl,
  title={{H}omotopy{C}ontinuation.jl: {A} {P}ackage for {H}omotopy {C}ontinuation in {J}ulia},
  author={Breiding, Paul and Timme, Sascha},
  booktitle={International Congress on Mathematical Software},
  pages={458--465},
  year={2018},
  organization={Springer}
}

A preprint of this paper is freely available.

homotopycontinuation.jl's Issues

Monodromy solver support and trace test

We should deliver an official monodromy implementation. This should work with either a parameter homotopy or an positive dimensional system. Additional an implementation of the trace test would be nice. This would enable us to verify irreducibility of varieties as well as computing the degree of a variety.

Error tagging new release

The tag name "v0.2.0-alpha.2" is not of the appropriate SemVer form (vX.Y.Z).
cc: @saschatimme

Add support for callbacks

There are often use cases where we look for one specific solution, e.g., the existence of a real solution with certain attributes. It would be nice that the user could supply a function which is called for each (finite) result such that the computation stops if the function returns true. This needs #124.

Better parameter homotopies / Symbolic tools

I think something along the lines of ParameterizedFunctions.jl would be really nice for parameter homotopies since currently no Julia packages has automatic differentiation for complex input.

Support for polyhedral homotopies

It would be amazing to support polyhedral homotopies. Before we can plan exactly the integration we need to have a mixed volume compuation available. The best algorithm for this is Anders Jensen's Tropical Homotopy Continuation. Also note his implementation details in An Implementation of Exact Mixed Volume Computation.

Replace FixedPolySystem and extract Homotopy

I create a similar package FixedPolynomials which I plan to publish soon. It features a brand new evaluation algorithm which is up to seven times faster in my benchmarks (and even more for polynomials with higher degrees).
I also think we should extract the Homotopy module in a standalone package under the JuliaHomotopyContinuation organization.

It is probably the best to make these two changes at once. My plan is to get to this within the next 1-2 weeks.

Public API for Newton's method

It would be nice to have a simple method available to call Newton's method with a certain amount of iterations, etc.

Linear intersection systems

We should have a custom homotopy type (and maybe also custom system type) which represents the intersection of a polynomial system (or homotopy) with a linear space. This can happen by adding addtional linear equations or computing the solution space (A, b) of the defining equations and then solving H(Ax+b).

Error tagging new release

The tag name "v0.2.0-alpha.0" is not of the appropriate SemVer form (vX.Y.Z).
cc: @saschatimme

Error tagging new release

The tag name "untagged-8413ad4e6fc4648ef972" is not of the appropriate SemVer form (vX.Y.Z).
cc: @saschatimme

differentiate exported but not defined

The guide makes use of a function differentiate and this is also exported by the package, however, the function is not defined anywhere.

Add coefficient parameter homotopy.

Now we have a fast parameter homotopy, but it would also be nice to special case the parameter homotopy which simply changes all coefficients of a polynomial.

Implement SCLGEN preprocessing

Alex Morgan developed a preprocessing heuristic to improve the numerical behaviour of polynomial systems. It is described in Section 5 of his book Solving Polynomial Systems Using Continuation for Engineering and Scientific Problems, here is a Google books link. This should be optional.

Add conversions and promotions

We currently have no conversions and promotions for the homotopies. One problem is that MPoly also doesn't have any. So I think the best would to wait until we replaced that.

Monodromy should check for solutions at infinity

It can happen (by path jumping) that we end up with a nonsingular solution at infinity. This completely kills the monodromy routine. We need a good check to verify that a solution is not at infinity.

Possible parallelism not working in Jupyter (or other slowdown as compared to Juno).

Been looking around some more. I get the impression that this uses some form of parallel computing? At least when is use "top" in ubuntu the %CPU regularly surpasses 100%.

However, this only seems to happen in Juno (in atom), but not in some Jupyter version (notebook or lab). I again have this sample code:

using HomotopyContinuation
@polyvar w w2 w2v v w2v2 vP σB w2σB vPp phos
poly = [(-1 * 0.7 * w + -2 * 3600.0 * (w ^ 2 / 2) + 2 * 18.0 * w2)*(0.2 + σB) + 4.0 * 0.4 * (1 + 30.0σB),
 -1 * 0.7 * w2 + 3600.0 * (w ^ 2 / 2) + -1 * 18.0 * w2 + -1 * 3600.0 * w2 * v + 18.0w2v + 36.0w2v + -1 * 3600.0 * w2 * σB + 18.0w2σB,
 -1 * 0.7 * w2v + 3600.0 * w2 * v + -1 * 18.0 * w2v + -1 * 3600.0 * w2v * v + 18.0w2v2 + -1 * 36.0 * w2v + 36.0w2v2 + 1800.0 * w2σB * v + -1 * 1800.0 * w2v * σB,
 (-1 * 0.7 * v + -1 * 3600.0 * w2 * v + 18.0w2v + -1 * 3600.0 * w2v * v + 18.0w2v2 + -1 * 1800.0 * w2σB * v + 1800.0 * w2v * σB + 180.0vPp)*(0.2 + σB) + 4.5 * 0.4 * (1 + 30.0σB),
 -1 * 0.7 * w2v2 + 3600.0 * w2v * v + -1 * 18.0 * w2v2 + -1 * 36.0 * w2v2,
 -1 * 0.7 * vP + 36.0w2v + 36.0w2v2 + -1 * 3600.0 * vP * phos + 18.0vPp,
 (-1 * 0.7 * σB + -1 * 3600.0 * w2 * σB + 18.0w2σB + 1800.0 * w2σB * v + -1 * 1800.0 * w2v * σB)*(0.2 + σB) + 0.4 * (1 + 30.0σB),
 -1 * 0.7 * w2σB + 3600.0 * w2 * σB + -1 * 18.0 * w2σB + -1 * 1800.0 * w2σB * v + 1800.0 * w2v * σB,
 -1 * 0.7 * vPp + 3600.0 * vP * phos + -1 * 18.0 * vPp + -1 * 180.0 * vPp,
 (phos + vPp) - 2.0]
@time res = HomotopyContinuation.solve(poly)

when I run this in Juno it takes about 6 seconds (was about 1 second a few weeks ago when I showed the previous problem, not sure when and why it changed). The top command show %CPU hitting 600 something. However, if I run the same code in JupyterLab the runtime is about 20 seconds and the CPU% never passes 100.

However, parallelism in general seems to work. I try this code

using Random
@time rand(6000,6000)*rand(6000,6000)

which runs in about 6 seconds in both Juno and Jupyter, and in both cases top shows that CPU% hits 600, well over 100.
(This uses HomotopyContinuation 0.3.2 and I should admit that I am no great expert in parallel computing)

Change order of computations in solve

Currently a run of solve consists of multiple stages. First we track all paths from 1 to the the beginning of the endgame zone, let's say 0.1.Then we make a path crossing check with possible reruns and then we play the endgame for each. The final stage is assembly of the output, refining if necessary etc.This leads to multiple problem. For high number of paths (in the millions) this leads to severe memory problems. Although CPU is no problem (this just takes a couple of hours) we run out of memory. But this is completely avoidable, since we could only hold a certain amount of paths in memory and store the rest on disk.To resolve this we should change the design such that each path is completely run independently from start to finish. An optional path crossing check could then be done at the end. This is also simplifies multi-node parallelism.

Replace MPoly

We should replace MPoly with something better. I think MultivariatePolynomials.jl is the best option, we just need to wait a little bit until there things settle down.

Start target homotopy with SVector as start solutions doesn't work

There is no PathResult constructor...

Error tagging new release

The tag name "v0.1" is not of the appropriate SemVer form (vX.Y.Z).
cc: @saschatimme

Store intermediate results automatically

For cluster support it would be nice to save intermediate results for long running computations automatically as well as give the possibilities to rerun an application.

Port alphaCertified to Julia

The program alphaCertified implements algorithms based on alpha-theory to certify solutions to polynomial systems using both exact rational arithmetic and arbitrary precision floating point arithmetic. It also implements an algorithm to certify whether a given point corresponds to a real solution to a real polynomial system, as well as algorithms to heuristically validate solutions to overdetermined systems.

Enable to solve a subset of all possible start solutions

Sometimes you just want to solve a small subset of all possible solutions to check certain characteristics of your system. The best way to enable is probably to publicly export the Solving.solver_startsolutions primitive and to export it properly.

Bug after upgrade to new FixedPolynomials version (master)

After updating FixedPolynomials I get the following error message:

using HomotopyContinuation, DynamicPolynomials

@polyvar x
solve([x-2])

MethodError: Cannot `convert` an object of type Array{FixedPolynomials.Polynomial{Complex{Float64}},1} to an object of type FixedPolynomials.JacobianConfig
This may have arisen from a call to the constructor FixedPolynomials.JacobianConfig(...),
since type constructors fall back to convert methods.
in solve at HomotopyContinuation/src/solve.jl:156
in #solve#82 at HomotopyContinuation/src/solve.jl:157
in #solve#92 at HomotopyContinuation/src/solve.jl:185
in #Solver#21 at HomotopyContinuation/src/solver.jl:109
in HomotopyContinuation.Pathtracker at HomotopyContinuation/src/pathtracker/type.jl:179
in #Pathtracker#5 at HomotopyContinuation/src/pathtracker/type.jl:179
in #Pathtracker#4 at HomotopyContinuation/src/pathtracker/type.jl:148
in Homotopies.PolynomialHomotopyConfig at Homotopies/src/config.jl:32

Collect more performance information

In order to better analyse the performance of our code it would be good to collect more metrics, in particular

• Number of functions evaluations
• Number of Jacobian evaluations
• Number of linear systems solved
• Number of accepted steps
• Number of rejected steps

How to track a path (for bifurcation analysis)

I am looking to track how a solution changes as we change a parameter of the system. The idea would be to solve the system for a certain set of parameters, then use that solution as a base to solve the same system, but with one parameter changed. Then I would extract that path between the two solutions and use that to create a bifurcation diagram. There are some info on path tracking in the documentation, but I have so far been unable to get it to work.

I have been attempting this:

using DynamicPolynomials, HomotopyContinuation
@polyvar x[1:2] c[1:2]
f = [x[1] - x[2] - c[1], c[2]*x[1]-x[2]-1]

p_start = [2.,2.]
f_start = subs.(f, Ref(c=> p_start))
result_start = solve(f_start, report_progress=false)

Now I can find a solution for new parameters with

p =  [2.,10.]
result = solve(f, solutions(result_start), parameters=c, p₁=p_start, p₀=p)

however, this does not seem to give me path info, so I try to use the pathtracking options. I can use:

path = pathtracker_startsolutions(f, solutions(result_start), parameters=c, p₁=p_start, p₀=p)

but it is here that the documentation starts to become hard to understand, exactly what do I do with this one. I try

HomotopyContinuation.PathTracking.PathTracker(path)
HomotopyContinuation.PathTracking.PathTracker(path[1].homotopy,result_path[2],1.,0.)
HomotopyContinuation.PathTracking.track!(path)
HomotopyContinuation.PathTracking.track!(path[1])
HomotopyContinuation.PathTracking.track!(path[1], path[2], 1., 0.)
HomotopyContinuation.PathTracking.track!(path[2],path[1], solutions(result), 1., 0.)
HomotopyContinuation.PathTracking.PathTrackerResult(path)

but all returns MethodError: no method matching (and if I do not start with HomotopyContinuation.PathTracking it does not work at all).

What odes work is

HomotopyContinuation.PathTracking.PathTrackerResult(path[1])

but it does not seem to give me a path, only a solution.

I get the impression that this is doable, and it might just be that I have been unable to understand the documentation page (https://www.juliahomotopycontinuation.org/HomotopyContinuation.jl/stable/pathtracking/#Path-tracking-1).

Some queries from a new users

I was happily surprised to discovered that homotopy continuation was implemented in Julia. I have been trying the package out, however there are a few non obvious things which I do not understand (this might partially be due to my lack of experience in the area). I am unsure whenever this is the right place to ask my questions, but it was the best I found (if not I apologise).

First, it seems that when there is a singularity present an additional 1.0 is added to the solution array, e.g.

@polyvar x y z w 
f = [x*y+z^2, y-3x*z+z^2+2, x^3+y^2+z-1];
sols = solutions(solve(f), only_real=true)
for i = 1:length(sols)
   println(real(sols[i].solution))
end

outputs

[3.62126e-46, 1.50418e-25, 1.0, 8.62497e-7]
[-1.00432e-24, -7.17904e-17, 1.0, -7.61496e-7]
[-1.40608, 2.36735, -1.82447]
[-76.8682, 674.107, -227.634]

where both the first and second solution contains a singularity. However, I did not find anything on what the actual meaning of the 1.0 and its position is.

Next I have a few things which all might be regarding my lack of knowledge in the area, I am no expert. However, when I run the part

sols = solutions(solve(f), only_real=true)
for i = 1:length(sols)
   println(real(sols[i].solution))
end

several times I do not all time get the previously mentioned solution, sometimes the solution [-76.8682, 674.107, -227.634] is not found (roughly half of the time). In addition this solution does not actually seem to be a good solution, indeed

sol = [-76.8682, 674.107, -227.634]
for i = 1:length(f)
   println(f[i](x=> sol[1] , y=>sol[2], z=>sol[3]),"   ") 
end

returns

-0.15374140000494663   
-0.10256040000604116   
-1.070417910619625   

is the solving method simply inaccurate for certain cases, not running long enough, or am I doing something wrong? The documentations says

The aim of this project is twofold: establishing a fast numerical polynomial solver in Julia and at the same time providing a highly customizable algorithmic environment for researchers for designing and performing individual experiments.

Then the documentations starts with a simple example of solving a polynomial system (what I want to do and also what I am trying here), and then continue with advanced stuff of defining your own homotopies and pathtracking. I am uncertain where the "fast numerical polynomial solver" (what I am interested in) ends and the "highly customizable algorithmic environment for researchers for designing and performing individual experiments" starts.

If the issue lies in my ability to use the methods and that they require some experience, do you maybe have a suggested reading for understanding more?

A further example where I struggle. For 4 variable systems I have been unable to find proper solutions for even simple linear systems like

@polyvar x y z w 
f = [w+x-4y, z+3-2x, 3w+x, 10-w+y];
sols = solutions(solve(f), only_real=true)
for i = 1:length(sols)
   println(real(sols[i].solution))
end

This gives

[6.66667, -20.0, -3.33333, -43.0]

However,

sol = [6.66667, -20.0, -3.33333, -43.0]
for i = 1:length(f)
   println(f[i](x=> sol[1] , y=>sol[2], z=>sol[3], w=>sol[4]),"   ") 
end

outputs

43.666669999999996   
-13.66667   
-122.33333   
33.0

which is distinctively non zero.

Support for multiprojective homotopies

We currently only support simple projective space. But multiprojective problems arise quite often in applications, so we should support this. This would need an change of the AffinePatches interfaces, ProjectiveVectors as well as the end game.

Naming consistency

Right now we are all over the place in our naming. Couple apis have multiple words phrased together as one word (e.g. startparameters) and sometimes we use snake case (e.g. report_progress). We should clean this up prior to 1.0.

Use Literate.jl for examples

We should consider using https://github.com/fredrikekre/Literate.jl for our examples. This would allow us to easily provide executable Jupyter notebooks with no additional overhead.

Initial Release Meta Issue

Everything I want to be done before the initial release:

• Support for adaptive precision #14 . For this I started to port the QD library to julia in the form of
  HigherPrecision.jl. Double-double arithmetic is already implemented. I have to see whether I have the time for quad-double arithmetic.
• Revise FixedPolynomials.jl and the Homotopy.jl API (as a consequence of former). I realised that the performance in FixedPolynomials for computing jacobians is really suboptimal. This should result in a nice performance boost.
• Add path-crossing checks
• We need to detect solutions at inifinity properly. I already added the condition of the jacobian and the magnitude of the homogenous coordinate to the result
• We need to get serious with documentation

Update: HigherPrecision.jl is now released with Doube-double precision.

Compute with equivalence classes in monodromy

In monodromy_solve we can declare group actions on the solutions, which is really handy. For this we apply for each new solution the group action and check whether the resulting solutions are new solutions. This works great, but I was wondering whether it wouldn't be smarter to only store equivalence classes of group actions. With this I mean, we store one solution and not the results of the group actions.
The logic after returning to the main node would then be (pseudo code)

if !iscontained(node, y)
    if any(y_i -> iscontained(node, y_i), group_action(y))
        return nothing
    end
    add!(node, y)
end

This would have the advantage that we only have to track the "true" solution count.

Parameter homotopy only works with polynomial input

Our parameter homotopy setup currently doesn't work with AbstractSystems. There's no reason why this shouldn't work and soon I will have a speficic application where I actually need this.

Add show functions everywhere and improve existing Juno integration

Currently not all defined types have a custom show(io, x) function defined.
This is in particular annoying for development and use in the REPL. Therefore, it would be nice to fix this. Also the existing show functions sometimes define show(io, ::MIME"text/plain", x) instead of show(io, x). This was necessary since the Juno integration didn't play nice with overleading the simple show function. Luckily this got fixed now and we should follow the behaviour detailed here.

Add benchmarks

I think it would be valuable to setup a proper benchmark enviroment. With this we could track performance regressions as well as it would be easy to compare different algorithms.

Encode size of system in types

Our polynomial systems have a fixed size, so we could start to encode this size in the relevant types .
For example, since the vectors with which we deal during pathtracking have always the same size, we can get extra performance with https://github.com/JuliaArrays/StaticArrays.jl.

BoundsError

With the following code:

using HomotopyContinuation
using DynamicPolynomials
@polyvar x
function _solve(f)
     F = PolySystem(f)
     G, startsolutions = totaldegree(F)
     H = StraightLineHomotopy(G, F)
     map(s -> s.solution, solve(H, startsolutions))
end
_solve([x-1, x-2])

I expect to get no solution as a result since x cannot be both 1 and 2 but I get the following error:

ERROR: BoundsError: attempt to access 1×2 Array{Int64,2} at index [2, 1]
Stacktrace:
 [1] (::HomotopyContinuation.HomConBase.##1#3{Complex{Float64},Array{Int64,1},Array{Complex{Float64},1},Int64})(::Int64) at /home/blegat/.julia/v0.6/HomotopyContinuation/src/HomConBase/utilities.jl:23
 [2] collect_to!(::Array{FixedPolySystem.Poly{Complex{Float64}},1}, ::Base.Generator{UnitRange{Int64},HomotopyContinuation.HomConBase.##1#3{Complex{Float64},Array{Int64,1},Array{Complex{Float64},1},Int64}}, ::Int64, ::Int64) at ./array.jl:474
 [3] _collect(::UnitRange{Int64}, ::Base.Generator{UnitRange{Int64},HomotopyContinuation.HomConBase.##1#3{Complex{Float64},Array{Int64,1},Array{Complex{Float64},1},Int64}}, ::Base.EltypeUnknown, ::Base.HasShape) at ./array.jl:455
 [4] total_degree_helper(::Array{Symbol,1}, ::Array{Int64,1}, ::Array{Complex{Float64},1}) at /home/blegat/.julia/v0.6/HomotopyContinuation/src/HomConBase/utilities.jl:19
 [5] #totaldegree#5(::Bool, ::Function, ::Type{Complex{Float64}}, ::Array{Symbol,1}, ::Array{Int64,1}) at /home/blegat/.julia/v0.6/HomotopyContinuation/src/HomConBase/utilities.jl:54
 [6] #totaldegree#6(::Array{Any,1}, ::Function, ::FixedPolySystem.PolySystem{Int64}) at /home/blegat/.julia/v0.6/HomotopyContinuation/src/HomConBase/utilities.jl:63
 [7] _solve(::Array{DynamicPolynomials.Polynomial{true,Int64},1}) at ./REPL[11]:3

Reevaluate the number of samples in the Cauchy endgame

I think it would be beneficial to have a more dynamic logic managing the number of samples used to predict a point in the Cauchy endgame. Currently it is just one fixed value, while I think it should linearly (?) increase with the number of predictions. With this we do not do too much work for the values where we actually do not need an endgame (which usually only need 2 prediction) and where
the algorithm is uncertain we increase the number of samples.

Support SemialgebraicSets.jl

It would be nice to support SemialgebraicSets.jl.
To add the support should be straight forward since we only need to define a HomotopyContinuationSolver <: AbstractAlgebraicSolve.

I currently prefer to do this after SemialgebraicSets.jl is released, since we plan to release a first version and therefore shouldn't introduce another unreleased dependency.

cc @blegat

New Newton method for ill-conditioned systems

For paths going to infinity or to singular solutions the Jacobian becomes more and more ill-conditioned.
This reduces the convergence speed of Newton's method and decreases the accuracy of the prediction methods.
Ultimately this leads to a failure of the path tracker.

To overcome this we need to change the way we solve the linear systems arising in Newton's method and the predictor methods for ill-conditioned systems.

In A Modified Newton Method for the Solution
of Ill-Conditioned Systems of Nonlinear Equations with Application to Multiple Shooting in such a way as to apply to the solution of ill-conditioned systems of nonlinear equations, i.e. systems having a "nearly singular" Jacobian at some iterate.
The core idea is to use the Moore-Penrose Pseudo-inverse instead of the normal inverse (which mathematically still exists). The core ideas are explained in Section 2 and in Section 4 implementation details are given. In particular we do not need to implement a SVD but rather only a QR decomposition with column pivoting (qr(A, Val(true)) in Julia).

The specific tasks to do are:

• Implement the proposed method to solve ill-conditioned systems
• Determine a criterion when to switch to this implementation during the path tracking
• Implement an automatic switching

Properly handle exceptions

An exception during a single path run shouldn't crash the whole program.

error when some elements of f are zero

Hi, I'm on Julia 0.6.2 and I've had the following issue on the latest release and master of HomotopyContinuation. I am trying to find the roots of multivariate polynomials as a subroutine of my code, but when some of the derivatives in the gradient f are zero it triggers an error with an unhelpful error message. It would be convenient if the package still worked in this case.

using MultivariatePolynomials
using DynamicPolynomials
using HomotopyContinuation

n = 4
@polyvar x[1:n]
c = [1, 1, 0, 0] # sometimes c will contain no zeros
p = sum(c[i]*x[i]^2 for i in 1:n)
# x1^2 + x2^2
f= [differentiate(p, x[i]) for i in 1:n]
# DynamicPolynomials.Polynomial{true,Int64}[2x1, 2x2, 0, 0]
HomotopyContinuation.solve(f)

ERROR: BoundsError: attempt to access 4×0 Array{Int64,2} at index [Base.Slice(Base.OneTo(4)), 1]
Stacktrace:
 [1] throw_boundserror(::Array{Int64,2}, ::Tuple{Base.Slice{Base.OneTo{Int64}},Int64}) at ./abstractarray.jl:434
 [2] checkbounds at ./abstractarray.jl:362 [inlined]
 [3] macro expansion at ./multidimensional.jl:494 [inlined]
 [4] _getindex at ./multidimensional.jl:491 [inlined]
 [5] getindex at ./abstractarray.jl:883 [inlined]
 [6] degree(::FixedPolynomials.Polynomial{Complex{Float64}}) at /home/coey/.julia/v0.6/FixedPolynomials/src/poly.jl:146
 [7] macro expansion at ./broadcast.jl:155 [inlined]
 [8] macro expansion at ./simdloop.jl:73 [inlined]
 [9] macro expansion at ./broadcast.jl:149 [inlined]
 [10] _broadcast! at ./broadcast.jl:141 [inlined]
 [11] broadcast_t at ./broadcast.jl:270 [inlined]
 [12] broadcast_c at ./broadcast.jl:316 [inlined]
 [13] broadcast at ./broadcast.jl:455 [inlined]
 [14] #totaldegree_startsystem#52(::Bool, ::Function, ::Array{FixedPolynomials.Polynomial{Complex{Float64}},1}) at /home/coey/.julia/v0.6/Homotopies/src/totaldegree.jl:99
 [15] #totaldegree#51(::Array{Any,1}, ::Function, ::Type{Homotopies.StraightLineHomotopy{Complex{Float64}}}, ::Array{FixedPolynomials.Polynomial{Complex{Float64}},1}) at /home/coey/.julia/v0.6/Homotopies/src/totaldegree.jl:75
 [16] #totaldegree#49(::Array{Any,1}, ::Function, ::Type{Homotopies.StraightLineHomotopy{Complex{Float64}}}, ::Array{DynamicPolynomials.Polynomial{true,Int64},1}) at /home/coey/.julia/v0.6/Homotopies/src/totaldegree.jl:67
 [17] #solve#82(::Array{Any,1}, ::Function, ::Array{DynamicPolynomials.Polynomial{true,Int64},1}) at /home/coey/.julia/v0.6/HomotopyContinuation/src/solve.jl:156
 [18] solve(::Array{DynamicPolynomials.Polynomial{true,Int64},1}) at /home/coey/.julia/v0.6/HomotopyContinuation/src/solve.jl:156

Adapt precision

Bertini uses higher precision numbers during the endgame and it would be nice if we would also do something like that.
To geht arbitrary precision we could use BigFloat, but that is super slow. Another possibility would be to switch to DoubleDouble, which would give at least double the precision.

The code should be written so that we would probably only need some conversions to integrate this.

Possible decrease in runtime since earlier version?

This may or may not be an actual issue. However, I have been trying to solve a system, which I previously remember solving a few months ago (think this was still on late stage 0.6, possibly around the arrival of 1.0). However, now when I solve the same system the runtime has increased a lot. While I remember it to be around 1 second, it is now around 150 seconds.

I am looking at the following system:

using HomotopyContinuation
@polyvar w w2 w2v v w2v2 vP σB w2σB vPp phos
poly = [(-1 * 0.7 * w + -2 * 3600.0 * (w ^ 2 / 2) + 2 * 18.0 * w2)*(0.2 + σB) + 4.0 * 0.4 * (1 + 30.0σB),
 -1 * 0.7 * w2 + 3600.0 * (w ^ 2 / 2) + -1 * 18.0 * w2 + -1 * 3600.0 * w2 * v + 18.0w2v + 36.0w2v + -1 * 3600.0 * w2 * σB + 18.0w2σB,
 -1 * 0.7 * w2v + 3600.0 * w2 * v + -1 * 18.0 * w2v + -1 * 3600.0 * w2v * v + 18.0w2v2 + -1 * 36.0 * w2v + 36.0w2v2 + 1800.0 * w2σB * v + -1 * 1800.0 * w2v * σB,
 (-1 * 0.7 * v + -1 * 3600.0 * w2 * v + 18.0w2v + -1 * 3600.0 * w2v * v + 18.0w2v2 + -1 * 1800.0 * w2σB * v + 1800.0 * w2v * σB + 180.0vPp)*(0.2 + σB) + 4.5 * 0.4 * (1 + 30.0σB),
 -1 * 0.7 * w2v2 + 3600.0 * w2v * v + -1 * 18.0 * w2v2 + -1 * 36.0 * w2v2,
 -1 * 0.7 * vP + 36.0w2v + 36.0w2v2 + -1 * 3600.0 * vP * phos + 18.0vPp,
 (-1 * 0.7 * σB + -1 * 3600.0 * w2 * σB + 18.0w2σB + 1800.0 * w2σB * v + -1 * 1800.0 * w2v * σB)*(0.2 + σB) + 0.4 * (1 + 30.0σB),
 -1 * 0.7 * w2σB + 3600.0 * w2 * σB + -1 * 18.0 * w2σB + -1 * 1800.0 * w2σB * v + 1800.0 * w2v * σB,
 -1 * 0.7 * vPp + 3600.0 * vP * phos + -1 * 18.0 * vPp + -1 * 180.0 * vPp,
 (phos + vPp) - 2.0]
@time res = HomotopyContinuation.solve(poly)

Admittedly I am not expert on this, but I figured I should raise the issue to those who know more. Does the runtime seem plausible? I might simply be misremembering things, but I think I remember right, and that if it indeed took ~150s when I did this earlier I would probably remember it.

Thanks

P.S.
(on a more general note, is there any, not to heavy, resource introducing Homotopy Continuation as a method to solve polynomial systems which you would recommend? Have background in mathematics.)

Vendor in DynamicPolynomials or TypedPolynomials

I wonder whether we should already include DynamicPolynomials or TypedPolynomials, since we use them in every example....

Large memory footprint when a lot of paths need to be tracked.

I currently run a calculation with total degree 6^9 (about 10 million) and the needed memory for the start solutions alone is 14gb which seems way too high. It seems that the Vector{Vector{Complex128}} representation has a lot of overhead. So maybe it would be useful to have the start solutions and the result stuff as a SVector.

system=SPSystem is a terrible API

Nobody besides me and Paul knows probably of this. We should keep this API but also add a new one like fast=true|false|nothing (default=nothing) where for nothing we decide things based on a heuristic.

Clean up unique points implementations

The new monodromy solver (#142) also introduces a new data structure (UniquePoints) for detecting duplicate points. This works much better than our existing solution used in the path crossing check. We should therefore throw out the old code and only us the new infrastructure.

Make pathcrossing check optional

The path crossing check can be rather expensive and can lead to memory problems. Therefore it should be possible to disable this.

Visualization of paths

A requested feature is the support of visualization of single paths. The idea is to have the complex plane and to plot each coordinate in an different color. This can be accomplished using the pathtracker_startsolutions function to assemble a pathtracker. The best it probably to use RecipesBase.jl to implement the plotting capabilities without directly depending on a plotting package.