Decoupled Homotopy

[oameye]

https://arxiv.org/pdf/[2208.08179v1](https://arxiv.org/pdf/2208.08179v1)

using HomotopyContinuation, LinearAlgebra
import Base.Iterators: product


struct OscillatorCache{T₁,T₂<:Tracker}
    s::Vector{ComplexF64}
    tracker_step_1::T₁
    tracker_step_2::T₂
end

## Define N = number of oscillators
N = 2
M = 2(N + 3)N
println("N = $N")
println("5^N = $(5^N)")

## Start Parameters p₀
p₀ = [randn(ComplexF64, 8) for _ in 1:N]

## Middle Parameters p₁
p₁ = [zeros(8 * N); randn(ComplexF64, M - N * 8)]

## Target Parameters p
p = randn(Float64, M)

## Start System
time_ss = @elapsed begin
    println("\n Computing Start System")
    ## Start System
    @var u₀, v₀
    @var a₀[1:4], b₀[1:4]
    n² = u₀^2 + v₀^2

    f = a₀[1] * n² * u₀ + a₀[2] * u₀ + a₀[3] * v₀ + a₀[4]
    g = b₀[1] * n² * v₀ + b₀[2] * u₀ + b₀[3] * v₀ + b₀[4]
    variables₀ = [u₀; v₀]
    parameters₀ = [a₀; b₀]

    Start_System = System([f; g],
        variables=variables₀,
        parameters=parameters₀)

    ## Start solutions
    p_single = randn(ComplexF64, 8)
    S₀ = solve(Start_System, target_parameters=p_single,
        show_progress=false)

    S = solve(Start_System,
        solutions(S₀),
        start_parameters=p_single,
        target_parameters=p₀,
        transform_result=(R, pᵢ) -> solutions(R),
        show_progress=false)

end
println("setting up start system = $time_ss seconds\n")

## Target system
@var u[1:N], v[1:N]
@var a[1:4, 1:N], b[1:4, 1:N]
@var c[1:N, 1:N], d[1:N, 1:N]
params = [vcat([[a[:, i]; b[:, i]] for i in 1:N]...);
    [c[i, j] for i in 1:N, j = 1:N if i != j];
    [d[i, j] for i in 1:N, j = 1:N if i != j]]

C = map(product(1:N, 1:N)) do (i, j)
    if i == j
        0
    else
        c[i, j]
    end
end
D = map(product(1:N, 1:N)) do (i, j)
    if i == j
        0
    else
        d[i, j]
    end
end

fs = [subs(f,
    variables₀ => [u[i], v[i]],
    parameters₀ => [a[:, i]; b[:, i]]) for i in 1:N]
gs = [subs(g,
    variables₀ => [u[i], v[i]],
    parameters₀ => [a[:, i]; b[:, i]]) for i in 1:N]

F = System([fs + C * v; gs + D * u],
    variables=[u; v],
    parameters=params)

## Tracking
@var t

F₁ = System(F([u; v], [vcat(p₀...); zeros(M - N * 8)] + (1 - t) .* p₁),
    variables=[u; v], parameters=[t])
F₂ = System(F([u; v], t .* p₁ + (1 - t) .* p),
    variables=[u; v], parameters=[t])
tracker₁ = Tracker(ParameterHomotopy(F₁, [1], [0]))
tracker₂ = Tracker(ParameterHomotopy(F₂, [1], [0]))

start_solution_iterator = Iterators.product(S...)

println("\n Solving Target System")
function oscillator_track(x::NTuple{N,Vector{ComplexF64}}, cache::OscillatorCache)
    s = cache.s
    for (j, xᵢ) in enumerate(x)
        s[j] = xᵢ[1]
        s[N+j] = xᵢ[2]
    end
    T = cache.tracker_step_1
    track!(T, s, 1, 0)
    s = T.state.x
    T = cache.tracker_step_2
    res = track!(T, s, 1, 0)
    if is_success(res)
        return copy(T.state.x)
    end
    nothing
end

cache = OscillatorCache(zeros(ComplexF64, 2 * N), tracker₁, tracker₂)

X = Vector{Union{Vector{ComplexF64},Nothing}}()
time_b = @elapsed for x in start_solution_iterator
    push!(X, oscillator_track(x, cache))
end
println("solving target system = $time_b seconds\n")
success_rate = 1 - count(X .== nothing) / 5^N
println("success rate = $success_rate\n")


#

using HomotopyContinuation, LinearAlgebra
import Base.Iterators: product

struct OscillatorCache{T<:Tracker}
    s::Vector{ComplexF64}
    tracker::T
end


## Define N = number of oscillators
N = 2
M = 2(N + 3)N
println("N = $N")
println("5^N = $(5^N)")

## Target Parameters p
p = randn(Float64, M)

## Start Parameters p₀
p₀ = [p[8*(i-1)+1:8*i] for i in 1:N]

## Start System
time_s = @elapsed begin
    println("\n Computing Start System")
    ## Start System
    @var u₀, v₀
    @var a₀[1:4], b₀[1:4]
    n² = u₀^2 + v₀^2

    f = a₀[1] * n² * u₀ + a₀[2] * u₀ + a₀[3] * v₀ + a₀[4]
    g = b₀[1] * n² * v₀ + b₀[2] * u₀ + b₀[3] * v₀ + b₀[4]
    variables₀ = [u₀; v₀]
    parameters₀ = [a₀; b₀]

    Start_System = System([f; g],
        variables=variables₀,
        parameters=parameters₀)

    ## Start solutions
    p_single = randn(ComplexF64, 8)
    S₀ = solve(Start_System, target_parameters=p_single,
        show_progress=false)

    S = solve(Start_System,
        solutions(S₀),
        start_parameters=p_single,
        target_parameters=p₀,
        transform_result=(R, pᵢ) -> solutions(R),
        show_progress=false)

end
println("setting up start system = $time_s seconds\n")

## Target system
@var u[1:N], v[1:N]
@var a[1:4, 1:N], b[1:4, 1:N]
@var c[1:N, 1:N], d[1:N, 1:N]
params = [vcat([[a[:, i]; b[:, i]] for i in 1:N]...);
    [c[i, j] for i in 1:N, j = 1:N if i != j];
    [d[i, j] for i in 1:N, j = 1:N if i != j]]

C = map(product(1:N, 1:N)) do (i, j)
    if i == j
        0
    else
        c[i, j]
    end
end
D = map(product(1:N, 1:N)) do (i, j)
    if i == j
        0
    else
        d[i, j]
    end
end

fs = [subs(f,
    variables₀ => [u[i], v[i]],
    parameters₀ => [a[:, i]; b[:, i]]) for i in 1:N]
gs = [subs(g,
    variables₀ => [u[i], v[i]],
    parameters₀ => [a[:, i]; b[:, i]]) for i in 1:N]

F = System([fs + C * v; gs + D * u],
    variables=[u; v],
    parameters=params)


## Tracking
@var t
q = p[8*N+1:M]
F = System(F([u; v], [vcat(p₀...); (1 - t) .* q]),
    variables=[u; v], parameters=[t])
tracker = Tracker(ParameterHomotopy(F, [1], [0]))

start_solution_iterator = Iterators.product(S...)

println("\n Solving Target System")

function oscillator_track(x::NTuple{N,Vector{ComplexF64}}, cache::OscillatorCache)
    s = cache.s
    for (j, xᵢ) in enumerate(x)
        s[j] = xᵢ[1]
        s[N+j] = xᵢ[2]
    end
    T = cache.tracker
    track!(T, s, 1, 0)
    res = track!(T, s, 1, 0)
    if is_success(res)
        return copy(T.state.x)
    end
    nothing
end

cache = OscillatorCache(zeros(ComplexF64, 2 * N), tracker)

X = Vector{Union{Vector{ComplexF64},Nothing}}()
time_b = @elapsed for x in start_solution_iterator
    push!(X, oscillator_track(x, cache))
end
println("solving target system = $time_b seconds\n")
success_rate = 1 - count(X .== nothing) / 5^N
println("success rate = $success_rate\n")