using HomotopyContinuation, DynamicPolynomials
@polyvar x y z PD p[1:3,1:3] T[1:2]

f1 = PD + p[1,1] + p[1,2]*T[1] + p[1,3] * T[1]^2 - x
# note that we eliminate the denominator in equation 2 and 3
f2 = (PD + p[2,1] + p[2,2]*T[1])*T[2] + p[2,3] - y * T[2]
f3 = (PD + p[3,1] + p[3,2]*T[1])*T[2] + p[3,3] - z * T[2]

f = [f1, f2, f3]

params = [vec(p); x; y; z]
generic_params = randn(ComplexF64, length(params)) # Note we should sample them randomly
f_generic = subs.(f, Ref(params => generic_params))

generic_result = solve(f_generic)

generic_solutions = solutions(generic_result)

# Construct a path tracker to track the two solutions to our specific instances
# p₁ and 
tracker = pathtracker(f, generic_solutions;
            parameters=params,
            # we just use som dummy parameters for now
            targetparameters=generic_params,
            startparameters=generic_params,
            affine=true)

specific_params = randn(length(params));

params

# Setup the new parameters
set_parameters!(tracker.homotopy, generic_params, specific_params);

# track to the two solutions
r1 = track(tracker, generic_solutions[1])
r2 = track(tracker, generic_solutions[2])

# Lets check that the tracking was successfull for both
r1.returncode == PathTrackerStatus.success && r2.returncode == PathTrackerStatus.success

r1.x

r2.x

if maximum(abs ∘ imag, r1.x) < 1e-8
    println("First solution: ", real.(r1.x))
end
if maximum(abs ∘ imag, r2.x) < 1e-8
    println("Second solution: ", real.(r2.x))
end