diff --git a/benchmarks/bio-chemical-rection-networks.jl b/benchmarks/bio-chemical-rection-networks.jl index b6892818..3589c8ed 100644 --- a/benchmarks/bio-chemical-rection-networks.jl +++ b/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) diff --git a/benchmarks/cyclooctane.jl b/benchmarks/cyclooctane.jl index 19ea1eb9..1847c18b 100644 --- a/benchmarks/cyclooctane.jl +++ b/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 diff --git a/benchmarks/run_system_benchmark.jl b/benchmarks/run_system_benchmark.jl index 6f61ab30..c5cc4b30 100644 --- a/benchmarks/run_system_benchmark.jl +++ b/benchmarks/run_system_benchmark.jl @@ -50,4 +50,4 @@ function run_benchmark(file_name) BenchmarkTools.save(file_name, res) return res -end \ No newline at end of file +end diff --git a/benchmarks/tritangents.jl b/benchmarks/tritangents.jl index 37ac5f97..7d481471 100644 --- a/benchmarks/tritangents.jl +++ b/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 diff --git a/src/solve.jl b/src/solve.jl index 23ebaac9..166fb5bb 100644 --- a/src/solve.jl +++ b/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, @@ -695,6 +697,7 @@ function solve( n = length(target_parameters) progress = show_progress ? make_many_progress(n; delay = 0.3) : nothing + many_solve( S, starts, @@ -749,10 +752,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 @@ -767,11 +777,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 diff --git a/test/solve_test.jl b/test/solve_test.jl index d190b2eb..9d9617e2 100644 --- a/test/solve_test.jl +++ b/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