Merge pull request #600 from oameye/main

fix: fix total degree with many parameters threaded (#599)

Author: PBrdng

benchmarks/bio-chemical-rection-networks.jl

@@ -6,82 +6,135 @@ using HomotopyContinuation, DynamicPolynomials, LinearAlgebra
 @polyvar x[1:3] p[1:10]
 
 N = 1000
-results1_direct = fill(0,8)
-results1_template = fill(0,8)
-results2_direct = fill(0,8)
-results2_template = fill(0,8)
-results3_direct = fill(0,8)
-results3_template = fill(0,8)
-results4_direct = fill(0,8)
-results4_template = fill(0,8)
+results1_direct = fill(0, 8)
+results1_template = fill(0, 8)
+results2_direct = fill(0, 8)
+results2_template = fill(0, 8)
+results3_direct = fill(0, 8)
+results3_template = fill(0, 8)
+results4_direct = fill(0, 8)
+results4_template = fill(0, 8)
 
 @time for i = 1:N
-    pol = [ -x[1]*x[3]*p[3] - x[1]*p[2] + p[1],
-            x[1]*x[2]*x[3]*p[8]*p[9] + x[1]*x[3]*p[7]*p[8]*p[9] - x[2]*x[3]*p[5]*p[6] - x[2]*p[5]*p[6]*p[10] - x[2]*x[3]*p[4] - x[2]*p[4]*p[10],
-            x[2] + x[3] - 1.0 ]
-    p_vals = [0.04,0.04,1.,1.,10.0,0.,0.04,35.,.1,.04]
+    pol = [
+        -x[1] * x[3] * p[3] - x[1] * p[2] + p[1],
+        x[1] * x[2] * x[3] * p[8] * p[9] + x[1] * x[3] * p[7] * p[8] * p[9] -
+        x[2] * x[3] * p[5] * p[6] - x[2] * p[5] * p[6] * p[10] - x[2] * x[3] * p[4] -
+        x[2] * p[4] * p[10],
+        x[2] + x[3] - 1.0,
+    ]
+    p_vals = [0.04, 0.04, 1.0, 1.0, 10.0, 0.0, 0.04, 35.0, 0.1, 0.04]
     pol_pars = DynamicPolynomials.subs.(pol, Ref(p => p_vals))
-    sol = solutions(solve(pol_pars, show_progress=false))
-    results1_direct[length(sol)+1] +=1
+    sol = solutions(solve(pol_pars, show_progress = false))
+    results1_direct[length(sol)+1] += 1
 
     p_template = randn(ComplexF64, 10)
     f_template = DynamicPolynomials.subs.(pol, Ref(p => p_template))
-    result_template = solutions(solve(f_template, show_progress=false))
-    sol_again = solutions(solve(pol, result_template, parameters=p, p₁=p_template,
-            p₀=ComplexF64.(p_vals), show_progress=false))
+    result_template = solutions(solve(f_template, show_progress = false))
+    sol_again = solutions(
+        solve(
+            pol,
+            result_template,
+            parameters = p,
+            p₁ = p_template,
+            p₀ = ComplexF64.(p_vals),
+            show_progress = false,
+        ),
+    )
     results1_template[length(sol_again)+1] += 1
 
-    pol2 = [ -x[1]^5*x[2]*p[5] + x[1]^4*x[3]*p[6]*p[8] - x[1]*x[2]*p[3]^4*p[5] + x[3]*p[3]^4*p[6]*p[8] - x[1]^5*p[7] + x[1]^4*x[3]*p[4] - x[1]*p[3]^4*p[7] + x[3]*p[3]^4*p[4] + x[1]^4*p[1] + x[1]^4*p[2] + p[1]*p[3]^4,
-             -x[1]^5*x[2]*p[5] - x[1]*x[2]*p[3]^4*p[5] - x[1]^4*x[2]*p[7] + x[1]^4*x[3]*p[4] - x[2]*p[3]^4*p[7] + x[3]*p[3]^4*p[4] + x[1]^4*p[1] + x[1]^4*p[2] + p[1]*p[3]^4,
-             x[1]*x[2]*p[5] - x[3]*p[6]*p[8] - x[3]*p[4] - x[3]*p[7] ]
-    p2_vals = [0.005, 0.1, 2.8,  10, 100, 0.1, 0.01, 0.]
+    pol2 = [
+        -x[1]^5 * x[2] * p[5] + x[1]^4 * x[3] * p[6] * p[8] - x[1] * x[2] * p[3]^4 * p[5] + x[3] * p[3]^4 * p[6] * p[8] - x[1]^5 * p[7] +
+        x[1]^4 * x[3] * p[4] - x[1] * p[3]^4 * p[7] +
+        x[3] * p[3]^4 * p[4] +
+        x[1]^4 * p[1] +
+        x[1]^4 * p[2] +
+        p[1] * p[3]^4,
+        -x[1]^5 * x[2] * p[5] - x[1] * x[2] * p[3]^4 * p[5] - x[1]^4 * x[2] * p[7] +
+        x[1]^4 * x[3] * p[4] - x[2] * p[3]^4 * p[7] +
+        x[3] * p[3]^4 * p[4] +
+        x[1]^4 * p[1] +
+        x[1]^4 * p[2] +
+        p[1] * p[3]^4,
+        x[1] * x[2] * p[5] - x[3] * p[6] * p[8] - x[3] * p[4] - x[3] * p[7],
+    ]
+    p2_vals = [0.005, 0.1, 2.8, 10, 100, 0.1, 0.01, 0.0]
     pol2_pars = DynamicPolynomials.subs.(pol2, Ref(p => p2_vals))
-    sol2 = solutions(solve(pol2_pars, show_progress=false))
-    results2_direct[length(sol2)+1] +=1
+    sol2 = solutions(solve(pol2_pars, show_progress = false))
+    results2_direct[length(sol2)+1] += 1
 
     p2_template = randn(ComplexF64, 8)
     f2_template = DynamicPolynomials.subs.(pol2, Ref(p => p2_template))
-    solve_templ2 = solve(f2_template, show_progress=false)
+    solve_templ2 = solve(f2_template, show_progress = false)
     result2_template = solutions(solve_templ2)
-    sol2_again = solutions(solve(pol2, result2_template, precision=PRECISION_ADAPTIVE, parameters=p[1:8], p₁=p2_template, p₀=ComplexF64.(p2_vals), show_progress=false))
+    sol2_again = solutions(
+        solve(
+            pol2,
+            result2_template,
+            precision = PRECISION_ADAPTIVE,
+            parameters = p[1:8],
+            p₁ = p2_template,
+            p₀ = ComplexF64.(p2_vals),
+            show_progress = false,
+        ),
+    )
 
-    results2_template[length(sol2_again)+1] +=1
+    results2_template[length(sol2_again)+1] += 1
 
     pol3 = pol2
-    p3_vals = [0.005, 0.1, 2.8,  10, 100, 0.1, 0.01, 1.] #Last parameter is 1 and not 0. Here there should be 3 real solutions instead of 1.
+    p3_vals = [0.005, 0.1, 2.8, 10, 100, 0.1, 0.01, 1.0] #Last parameter is 1 and not 0. Here there should be 3 real solutions instead of 1.
     pol3_pars = DynamicPolynomials.subs.(pol3, Ref(p => p3_vals))
-    sol3 = solutions(solve(pol3_pars, precision=PRECISION_ADAPTIVE, show_progress=false))
-    results3_direct[length(sol3)+1] +=1
+    sol3 =
+        solutions(solve(pol3_pars, precision = PRECISION_ADAPTIVE, show_progress = false))
+    results3_direct[length(sol3)+1] += 1
 
     p3_template = randn(ComplexF64, 8)
     f3_template = DynamicPolynomials.subs.(pol3, Ref(p => p3_template))
-    result3_template = solutions(solve(f3_template, show_progress=false))
-    sol3_again = solutions(solve(pol3, result3_template, parameters=p[1:8], p₁=p3_template, p₀=ComplexF64.(p3_vals), show_progress=false))
-    results3_template[length(sol3_again)+1] +=1
+    result3_template = solutions(solve(f3_template, show_progress = false))
+    sol3_again = solutions(
+        solve(
+            pol3,
+            result3_template,
+            parameters = p[1:8],
+            p₁ = p3_template,
+            p₀ = ComplexF64.(p3_vals),
+            show_progress = false,
+        ),
+    )
+    results3_template[length(sol3_again)+1] += 1
 
-    pol4 = [-x[1]^3*p[3] - x[1]*p[2]^2*p[3] + x[1]^2*p[1]]
-    p4_vals = [1., 0.2, 1.];
+    pol4 = [-x[1]^3 * p[3] - x[1] * p[2]^2 * p[3] + x[1]^2 * p[1]]
+    p4_vals = [1.0, 0.2, 1.0]
     pol4_pars = DynamicPolynomials.subs.(pol4, Ref(p => p4_vals))
-    sol4 = solutions(solve(pol4_pars, show_progress=false))
-    results4_direct[length(sol4)+1] +=1
+    sol4 = solutions(solve(pol4_pars, show_progress = false))
+    results4_direct[length(sol4)+1] += 1
 
     p4_template = randn(ComplexF64, 3)
     f4_template = DynamicPolynomials.subs.(pol4, Ref(p => p4_template))
-    result4_template = solutions(solve(f4_template, show_progress=false))
-    sol4_again = solutions(solve(pol4, result4_template, parameters=p[1:3], p₁=p4_template, p₀=ComplexF64.(p4_vals), show_progress=false))
-    results4_template[length(sol4_again)+1] +=1
+    result4_template = solutions(solve(f4_template, show_progress = false))
+    sol4_again = solutions(
+        solve(
+            pol4,
+            result4_template,
+            parameters = p[1:3],
+            p₁ = p4_template,
+            p₀ = ComplexF64.(p4_vals),
+            show_progress = false,
+        ),
+    )
+    results4_template[length(sol4_again)+1] += 1
 end
 
 # These numbers should all coincide
 println("Results for the first system")
 println(results1_direct)
-println(results1_template,"\n")
+println(results1_template, "\n")
 println("Results for the second system")
 println(results2_direct)
-println(results2_template,"\n")
+println(results2_template, "\n")
 println("Results for the third system")
 println(results3_direct)
-println(results3_template,"\n")
+println(results3_template, "\n")
 println("Results for the fourth system")
 println(results4_direct)
 println(results4_template)

benchmarks/cyclooctane.jl

@@ -4,13 +4,13 @@ using HomotopyContinuation, LinearAlgebra, DynamicPolynomials
 c² = 2
 @polyvar z[1:3, 1:6]
 z_vec = vec(z)[1:17] # the 17 variables in a vector
-Z = [zeros(3) z[:,1:5] [z[1,6]; z[2,6]; 0] [√c²; 0; 0]] # the eight points in a matrix
+Z = [zeros(3) z[:, 1:5] [z[1, 6]; z[2, 6]; 0] [√c²; 0; 0]] # the eight points in a matrix
 
 # define the functions for cyclooctane:
-F1 = [(Z[:, i] - Z[:, i+1]) ⋅ (Z[:, i] - Z[:, i+1]) - c² for i in 1:7]
-F2 = [(Z[:, i] - Z[:, i+2]) ⋅ (Z[:, i] - Z[:, i+2]) - 8c²/3 for i in 1:6]
-F3 = (Z[:, 7] - Z[:, 1]) ⋅ (Z[:, 7] - Z[:, 1]) - 8c²/3
-F4 = (Z[:, 8] - Z[:, 2]) ⋅ (Z[:, 8] - Z[:, 2]) - 8c²/3
+F1 = [(Z[:, i] - Z[:, i+1]) ⋅ (Z[:, i] - Z[:, i+1]) - c² for i = 1:7]
+F2 = [(Z[:, i] - Z[:, i+2]) ⋅ (Z[:, i] - Z[:, i+2]) - 8c² / 3 for i = 1:6]
+F3 = (Z[:, 7] - Z[:, 1]) ⋅ (Z[:, 7] - Z[:, 1]) - 8c² / 3
+F4 = (Z[:, 8] - Z[:, 2]) ⋅ (Z[:, 8] - Z[:, 2]) - 8c² / 3
 f = [F1; F2; F3; F4]
 
 n = 2 # dimension of the cyclooctane variety

benchmarks/run_system_benchmark.jl

@@ -50,4 +50,4 @@ function run_benchmark(file_name)
 
     BenchmarkTools.save(file_name, res)
     return res
-end
+end

benchmarks/tritangents.jl

@@ -5,5 +5,5 @@ f = equations(tritangents())
 res = solve(f)
 println("# nonsingular: ", nnonsingular(res)) # should be 720
 
-res = solve(f; start_system=:polyhedral)
+res = solve(f; start_system = :polyhedral)
 println("# nonsingular: ", nnonsingular(res)) # should be 720

src/solve.jl

@@ -130,6 +130,7 @@ function solver_startsolutions(
 )
     !isnothing(seed) && Random.seed!(seed)
 
+
     used_start_system = nothing
     if start_subspace !== nothing
         if target_parameters !== nothing
@@ -471,6 +472,7 @@ function solve(
     else
         solver, starts = solver_startsolutions(args...; kwargs...)
     end
+
     if many_parameters
         solve(
             solver,
@@ -703,6 +705,7 @@ function solve(
     n = length(target_parameters)
 
     progress = show_progress ? make_many_progress(n; delay = 0.3) : nothing
+
     many_solve(
         S,
         starts,
@@ -757,10 +760,17 @@ function many_solve(
     threading::Bool,
     catch_interrupt::Bool,
 ) where {flatten}
+
+
+    if isa(starts, TotalDegreeStartSolutionsIterator)
+        @warn "Solving for many parameters with total degree start system is not recommended. Instead, one should use a two-step approach: first solve a system with generic parameters, and track its solutions to the desired parameters. See https://www.juliahomotopycontinuation.org/guides/many-systems/."
+    elseif isa(starts, PolyhedralStartSolutionsIterator)
+        @error "Solving for many parameters with polyhedral start system is not implemented and also not recommended. Instead, one should use a two-step approach: solve a system with generic parameters, and track its solutions to the desired parameters. See https://www.juliahomotopycontinuation.org/guides/many-systems/."
+    end
     q = first(many_target_parameters)
     target_parameters!(solver, transform_parameters(q))
     if threading
-        res = threaded_solve(solver, starts; catch_interrupt = false)
+        res = threaded_solve(solver, collect(starts); catch_interrupt = false)
     else
         res = serial_solve(solver, starts; catch_interrupt = false)
     end
@@ -775,11 +785,14 @@ function many_solve(
     k = 1
     m = length(starts)
     update_many_progress!(progress, results, k, m; flatten = flatten)
+
+
+
     try
         for q in Iterators.drop(many_target_parameters, 1)
             target_parameters!(solver, transform_parameters(q))
             if threading
-                res = threaded_solve(solver, starts; catch_interrupt = false)
+                res = threaded_solve(solver, collect(starts); catch_interrupt = false)
             else
                 res = serial_solve(solver, starts; catch_interrupt = false)
             end

test/solve_test.jl

@@ -583,5 +583,45 @@
             intrinsic = false,
         )
         @test all(r -> nsolutions(first(r)) == 2, result1)
+
+        @testset "Many parameters threaded" begin
+            @var u1, v1, ω, α, γ, λ, ω0
+
+            eqs = [
+                -u1 * ω^2 +
+                u1 * ω0^2 +
+                (3 / 4) * u1^3 * α +
+                (3 / 4) * u1 * v1^2 * α +
+                (-1 / 2) * u1 * λ * ω0^2 +
+                v1 * γ * ω,
+                -v1 * ω^2 + v1 * ω0^2 + (3 / 4) * v1^3 * α - u1 * γ * ω +
+                (3 / 4) * u1^2 * v1 * α +
+                (1 / 2) * v1 * λ * ω0^2,
+            ]
+
+            F = System(eqs, parameters = [ω, α, γ, λ, ω0], variables = [u1, v1])
+
+            input_array = [
+                [0.9, 1.0, 0.01, 0.01, 1.1],
+                [0.9105263157894737, 1.0, 0.01, 0.01, 1.1],
+                [0.9210526315789473, 1.0, 0.01, 0.01, 1.1],
+                [0.9315789473684211, 1.0, 0.01, 0.01, 1.1],
+                [0.9421052631578948, 1.0, 0.01, 0.01, 1.1],
+                [0.9526315789473684, 1.0, 0.01, 0.01, 1.1],
+            ]
+
+            generic_parameters = randn(ComplexF64, 5)
+
+            R0 = solve(F; target_parameters = generic_parameters, threading = true)
+            R1 = solve(
+                F,
+                solutions(R0);
+                start_parameters = generic_parameters,
+                target_parameters = input_array,
+                threading = true,
+            )
+
+            @test length(R1) == 6
+        end
     end
 end