# Generate hill climbing data¶

#=
using Pkg
=#
using HomotopyContinuation
using JuMP
using GLPK
using LinearAlgebra
using Symbolics
using Convex
using SCS
using DataFrames
using DelimitedFiles
using CSV
@var s, t, r # Circle centered at (s, t) with radius r
@var u[1:3], v[1:3] # Three points of tangency (u[1], v[1]), ..., (u[3], v[3])
@var a[1:6], b[1:6], c[1:6] # Three fixed conics with coefficients defined by a[1], ..., a[6]; ... ; c[1], ..., c[6]
(HomotopyContinuation.ModelKit.Variable[a₁, a₂, a₃, a₄, a₅, a₆], HomotopyContinuation.ModelKit.Variable[b₁, b₂, b₃, b₄, b₅, b₆], HomotopyContinuation.ModelKit.Variable[c₁, c₂, c₃, c₄, c₅, c₆])
f_1 = (u[1] - s)^2 + (v[1] - t)^2 - r # the point (u[1], v[1]) lies on the circle

f_2 = (u[2] - s)^2 + (v[2] - t)^2 - r # the point (u[2], v[2]) lies on the circle

f_3 = (u[3] - s)^2 + (v[3] - t)^2 - r # the point (u[3], v[3]) lies on the circle

f_4 = a[1]*u[1]^2 + a[2]*u[1]*v[1] + a[3]*v[1]^2 + a[4]*u[1] + a[5]*v[1] + a[6] # the point (u[1], v[1]) lies on the conic defined by coefficients a[1:6]

f_5 = b[1]*u[2]^2 + b[2]*u[2]*v[2] + b[3]*v[2]^2 + b[4]*u[2] + b[5]*v[2] + b[6] # the point (u[2], v[2]) lies on the conic defined by coefficients b[1:6]

f_6 = c[1]*u[3]^2 + c[2]*u[3]*v[3] + c[3]*v[3]^2 + c[4]*u[3] + c[5]*v[3] + c[6] # the point (u[3], v[3]) lies on the conic defined by coefficients c[1:6]

f_7 = det([differentiate(f_1, [u[1], v[1]]) differentiate(f_4, [u[1], v[1]])]) # Circle and conic are tangent at (u[1], v[1])

f_8 = det([differentiate(f_2, [u[2], v[2]]) differentiate(f_5, [u[2], v[2]])]) # Circle and conic are tangent at (u[2], v[2])

f_9 = det([differentiate(f_3, [u[3], v[3]]) differentiate(f_6, [u[3], v[3]])]) # Circle and conic are tangent at (u[3], v[3])

# Polynomial system in variables

paramVec = collect(Iterators.flatten([a,b,c]))

F = System([f_1, f_2, f_3, f_4, f_5, f_6, f_7, f_8, f_9], variables = [u[1], v[1], u[2], v[2], u[3], v[3], s, t, r], parameters = paramVec);

# Polynomial system in coefficients

varVec = [u[1], v[1], u[2], v[2], u[3], v[3], s, t, r]

G = System([f_1, f_2, f_3, f_4, f_5, f_6, f_7, f_8, f_9], variables = [a[1], a[2], a[3], a[4], a[5], a[6], b[1], b[2], b[3], b[4], b[5], b[6], c[1], c[2], c[3], c[4], c[5], c[6]], parameters = varVec)
System of length 9
18 variables: a₁, a₂, a₃, a₄, a₅, a₆, b₁, b₂, b₃, b₄, b₅, b₆, c₁, c₂, c₃, c₄, c₅, c₆
9 parameters: u₁, v₁, u₂, v₂, u₃, v₃, s, t, r

-r + (-s + u₁)^2 + (-t + v₁)^2
-r + (-s + u₂)^2 + (-t + v₂)^2
-r + (-s + u₃)^2 + (-t + v₃)^2
a₆ + u₁*a₄ + u₁^2*a₁ + v₁*a₅ + v₁^2*a₃ + u₁*v₁*a₂
b₆ + u₂*b₄ + u₂^2*b₁ + v₂*b₅ + v₂^2*b₃ + u₂*v₂*b₂
c₆ + u₃*c₄ + u₃^2*c₁ + v₃*c₅ + v₃^2*c₃ + u₃*v₃*c₂
-2*(a₄ + 2*u₁*a₁ + v₁*a₂)*(-t + v₁) + 2*(a₅ + u₁*a₂ + 2*v₁*a₃)*(-s + u₁)
-2*(b₄ + 2*u₂*b₁ + v₂*b₂)*(-t + v₂) + 2*(b₅ + u₂*b₂ + 2*v₂*b₃)*(-s + u₂)
-2*(c₄ + 2*u₃*c₁ + v₃*c₂)*(-t + v₃) + 2*(c₅ + u₃*c₂ + 2*v₃*c₃)*(-s + u₃)
# Jacobian of polynomial system with respect to variables

function jacobian_vars(a, u)
return jacobian(F, u, a)
end

# Jacobian of polynomial system with respect to parameters

function jacobian_params(a, u)
return jacobian(G, a, u)
end

# Block Jacobian -- separated real and imaginary part for polynomial system w.r.t. variables

function block_jacobian_vars(a, u)
J = jacobian(F, u, a)
M = zeros(18, 18)
M[1:9, 1:9] = real(J)
M[10:18, 10:18] = imag(J)
return M
end

function block_jacobian_params(a, u)
J = jacobian(G, a, u)
M = zeros(18, 18)
M[1:9, 1:18] = real(J)
M[10:18, 1:18] = imag(J)
return M
end
block_jacobian_params (generic function with 1 method)
function complex_roots(root, a, real_sols, complex_sols)
p_max = 10^-1

lr = length(r_sols)
lc = length(c_sols)

model = Model(GLPK.Optimizer)

@variable(model, -p_max <= d_a[1:18] <= p_max)
@variable(model, d_img[1:9])
@variable(model, d_real[1:lr, 1:9])
@variable(model, s[1:9] >= 0)
@variable(model, t[1:9] >= 0)
@variable(model, d_x[1:lc, 1:18])

@objective(model, Min, dot(s+t, ones(9)))

@constraint(model, s - t .== d_img)

@constraint(model, imag(jacobian_vars(a, root))*d_img + imag(jacobian_params(a, root))*d_a .== 0)

@constraint(model, c[i=1:lr, j=i:lr], dot(real_sols[i]-real_sols[j], d_real[i:i,1:end]-d_real[j:j, 1:end]) >= 0)

@constraint(model, d[i=1:lr], (jacobian_vars(a, real_sols[i])*d_real[i:i, 1:end]')[1:end] + jacobian_params(a, real_sols[i])*d_a .== 0)

#@constraint(model, e[i=1,lc], (block_jacobian_vars(a, root) * (d_x[i:i, 1:end])')[1:end] + block_jacobian_params(a, root)*d_a .== 0)

#@constraint(model, f[i=1,lc], collect(Iterators.flatten([real(complex_sols[i]), imag(complex_sols[i])])) * d_x[i:i, 1:end] .>= 0)

optimize!(model)

n_p = value.(d_a)+a
n_root = solutions(solve(F, root, start_parameters = a, target_parameters = n_p, tracker_options = TrackerOptions(automatic_differentiation = 8, max_steps = 1000000, parameters = :conservative, extended_precision = true), show_progress = false))
n_r_sols = real_solutions(solve(F, r_sols, start_parameters = a, target_parameters = n_p, tracker_options = TrackerOptions(automatic_differentiation = 8, max_steps = 1000000, parameters = :conservative, extended_precision = true), show_progress = false))
n_c_sols = solutions(solve(F, c_sols, start_parameters = a, target_parameters = n_p, tracker_options = TrackerOptions(automatic_differentiation = 8, max_steps = 1000000, parameters = :conservative, extended_precision = true), show_progress = false))

return n_root[1], n_p, n_r_sols, n_c_sols
end
complex_roots (generic function with 1 method)
function real_roots(root1, root2, a, real_sols, complex_sols)

lr = length(real_sols)
lc = length(complex_sols)

p_max = 10^-2

model = Model(GLPK.Optimizer)

@variable(model, -p_max <= d_a[1:18] <= p_max)
@variable(model, d_real[1:lr, 1:9])
@variable(model, d_x[1:lc, 1:18])

eqn = dot(root1 - root2, d_real[end-1:end-1, 1:end] - d_real[end:end, 1:end])

@objective(model, Max, eqn)

@constraint(model, c[i=1:lr, j=i:lr], dot(real_sols[i]-real_sols[j], d_real[i:i,1:end]-d_real[j:j, 1:end]) >= 0)

@constraint(model, d[i=1:lr], (jacobian_vars(a, real_sols[i])*d_real[i:i, 1:end]')[1:end] + jacobian_params(a, real_sols[i])*d_a .== 0) #constraint on the real roots already there

#@constraint(model, (jacobian_vars(a, root1)*d_real[end-1:end-1, 1:end]')[1:end] + jacobian_params(a, root1)*d_a .== 0) # constraint on the new real root "root1"

#@constraint(model, (jacobian_vars(a, root2)*d_real[end:end, 1:end]')[1:end] + jacobian_params(a, root2)*d_a .== 0) # constraint on the new real root "root2"

#@constraint(model, e[i=1,lc], (block_jacobian_vars(a, complex_sols[i]) * (d_x[i:i, 1:end])')[1:end] + block_jacobian_params(a, complex_sols[i])*d_a .== 0)

#@constraint(model, f[i=1,lc], collect(Iterators.flatten([real(complex_sols[i]), imag(complex_sols[i])])) * d_x[i:i, 1:end] .>= 0)

optimize!(model)

return value.(d_a)+a
end
real_roots (generic function with 1 method)
function one_loop(p)

S = solve(F, target_parameters = p, show_progress = false);

rindex = findall(s -> is_real(s) && is_nonsingular(s), S)
cindex = findall(s -> !is_real(s) && is_nonsingular(s), S)

r_sols = solutions(S[rindex])
r_sols = real.(r_sols)
c_sols = solutions(S[cindex])

root = c_sols[1]

while norm(imag(root), Inf) > 10^-8
try
global root, p, r_sols, c_sols = complex_roots(root, p, r_sols, c_sols)
catch
r1 = real(root) .+ 10^-5
r2 = real(root) .- 10^-5

q = Convex.Variable(18)

problem = minimize(sumsquares(p-q))

problem.constraints += q[1]*r1[1]^2 + q[2]*r1[1]*r1[2] + q[3]*r1[2]^2 + q[4]*r1[1] + q[5]*r1[2] + q[6] == 0
problem.constraints += q[1]*r2[1]^2 + q[2]*r2[1]*r2[2] + q[3]*r2[2]^2 + q[4]*r2[1] + q[5]*r2[2] + q[6] == 0
problem.constraints += q[7]*r1[3]^2 + q[8]*r1[3]*r1[4] + q[9]*r1[4]^2 + q[10]*r1[3] + q[11]*r1[4] + q[12] == 0
problem.constraints += q[7]*r2[3]^2 + q[8]*r2[3]*r2[4] + q[9]*r2[4]^2 + q[10]*r2[3] + q[11]*r2[4] + q[12] == 0
problem.constraints += q[13]*r1[5]^2 + q[14]*r1[5]*r1[6] + q[15]*r1[6]^2 + q[16]*r1[5] + q[17]*r1[6] + q[18] == 0
problem.constraints += q[13]*r2[5]^2 + q[14]*r2[5]*r2[6] + q[15]*r2[6]^2 + q[16]*r2[5] + q[17]*r2[6] + q[18] == 0
problem.constraints += (q[5] + r1[1]*q[2] + 2*r1[2]*q[3])*(-r1[7] + r1[1]) - (q[4] + 2*r1[1]*q[1] + r1[2]*q[2])*(-r1[8] + r1[2]) == 0
problem.constraints += (q[5] + r2[1]*q[2] + 2*r2[2]*q[3])*(-r2[7] + r2[1]) - (q[4] + 2*r2[1]*q[1] + r2[2]*q[2])*(-r2[8] + r2[2]) == 0
problem.constraints += (q[11] + r1[3]*q[8] + 2*r1[4]*q[9])*(-r1[7] + r1[3]) - (q[10] + 2*r1[3]*q[7] + r1[4]*q[8])*(-r1[8] + r1[4]) == 0
problem.constraints += (q[11] + r2[3]*q[8] + 2*r2[4]*q[9])*(-r2[7] + r2[3]) - (q[10] + 2*r2[3]*q[7] + r2[4]*q[8])*(-r2[8] + r2[4]) == 0
problem.constraints += (q[17] + r1[5]*q[14] + 2*r1[6]*q[15])*(-r1[7] + r1[5]) - (q[16] + 2*r1[5]*q[13] + r1[6]*q[14])*(-r1[8] + r1[6]) == 0
problem.constraints += (q[17] + r2[5]*q[14] + 2*r2[6]*q[15])*(-r2[7] + r2[5]) - (q[16] + 2*r2[5]*q[13] + r2[6]*q[14])*(-r2[8] + r2[6]) == 0

Convex.solve!(problem, SCS.Optimizer)

p_ = round.(Convex.evaluate(q), digits = 14)

S_ = solve(F, target_parameters = p_);

rindex_ = findall(t -> is_real(t) && is_nonsingular(t), S_)
cindex_ = findall(t -> !is_real(t) && is_nonsingular(t), S_)

r_sols_ = solutions(S[rindex_])
r_sols_ = real.(r_sols_)
c_sols_ = solutions(S[cindex_])

return real_roots(r1, r2, p_, r_sols_, c_sols_)
#return p_
end
end
end
one_loop (generic function with 1 method)
p = rand(18)
18-element Vector{Float64}:
0.8350343218551497
0.9315554638276398
0.1614064506258529
0.4694662068562462
0.7581078380254038
0.48472578290889423
0.967919675103355
0.439078144594901
0.2890125754252907
0.799002076987044
0.7141935277885828
0.08733424926183997
0.8713339470821817
0.923713377091316
0.7418333529584032
0.8223292289031411
0.38570348414840405
0.4290246345110599
one_loop(p)
[32mTracking 256 paths... 100%|█████████████████████████████| Time: 0:00:02[39m
[34m  # paths tracked:                  256[39m
[34m  # non-singular solutions (real):  184 (22)[39m
[34m  # singular endpoints (real):      0 (0)[39m
[34m  # total solutions (real):         184 (22)[39m
18-element Vector{Float64}:
0.8270733456439268
0.94219614418061
0.169922746782527
0.45949984064525
0.74827454887878
0.49477380618
-0.09610882316597
-0.24474507198574208
0.06563005853736525
0.40003970692362
0.49952603549290214
-0.04480402332297458
-0.16058462211299987
-0.004036538157498566
-0.08634859817849999
0.4647471979690812
0.07857668021063
0.3048641063231938
i=1
df = DataFrame(A=Vector[], B=Vector[])
while i <= 50000
try
p = rand(18)
push!(df, [p, one_loop(p)])
global i += 1
catch
global i += 0
end
end
CSV.write("hc_data8.csv", df)
[32mTracking 256 paths... 100%|█████████████████████████████| Time: 0:00:00[39m
[34m  # paths tracked:                  256[39m
[34m  # non-singular solutions (real):  184 (8)[39m
[34m  # singular endpoints (real):      0 (0)[39m
[34m  # total solutions (real):         184 (8)[39m