Generate hill climbing data

#=
using Pkg
Pkg.add("HomotopyContinuation")
Pkg.add("JuMP")
Pkg.add("GLPK")
Pkg.add("LinearAlgebra")
Pkg.add("Symbolics")
Pkg.add("Convex")
Pkg.add("SCS")
Pkg.add("DataFrames")
Pkg.add("DelimitedFiles")
Pkg.add("CSV")
=#
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}:
8350343218551497
9315554638276398
1614064506258529
4694662068562462
7581078380254038
48472578290889423
967919675103355
439078144594901
2890125754252907
799002076987044
7141935277885828
08733424926183997
8713339470821817
923713377091316
7418333529584032
8223292289031411
38570348414840405
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