using HomotopyContinuation, DynamicPolynomials, LinearAlgebra, StaticPolynomials

# f : Sym^2(C^d)^n  -> Sym^d(C^n)
# X = im(f)
d = 2
n = 4
M = binomial(d - 1 + 2, 2) # = dim( Sym^2(C^d) )
N = binomial(n - 1 + d, d) # = dim( Sym^d(C^n) )
D = n * M # = dim( Sym^2(C^d)^n )

@polyvar x[1:n]
@polyvar a[1:D]

# Construct symmetric matrices with variable entries
A = map(0:n-1) do ℓ
   Aᵢ  = []
   for i in 1:d
            k = ℓ * M + sum(d-j for j in 0:i-1)
            push!(Aᵢ, [zeros(i-1); a[k-(d-i):k]])
   end
   Aᵢ = hcat(Aᵢ...)
   (Aᵢ + transpose(Aᵢ)) ./ 2
end

# Compute all monomials of a given degree d
function segre(x, d)
    if d == 1
        return x
    else
        return kron(x,segre(x,d-1))
    end
end
# ν = Veronese map
ν = unique(segre(x,d))

# Compute explicit the map f
g = det(sum(x[i] .* A[i] for i in 1:n))
f = [coefficient(g, v, x) for v in ν]
f = [subs(fᵢ, a[D] => 1.0, x => zeros(n)) for fᵢ in f] # make it affine

# Compute the codimension of X = im(f)
J = differentiate(f,a[1:D-1])
a₀ = randn(ComplexF64, D-1)
J₀ = [j(a[1:D-1] => a₀) for j in J]
dimX = rank(J₀) # = dim( X )
C = N - dimX
dim_fibers = D - 1 - dimX

# Compute a point f₀ on X and a linear space {Q f = q} of codimension dimX through it
f_SP = StaticPolynomials.PolynomialSystem(f)
f₀ = f_SP(a₀)
Q = randn(ComplexF64, dimX, N)
q = Q * f₀
# Compute a linear space {R a = r} of codimension dim_fibers through a₀
R = randn(ComplexF64, dim_fibers, D-1)
r = R * a₀

# Computing the fiber
# Monodromy by moving the linear space
@polyvar K[1:D]
fiber_cap_K = [
    f - f₀;
    transpose(K[1:D-1]) * a[1:D-1] - K[D];
    R[2:dim_fibers, :] * a[1:D-1] - r[2:dim_fibers];
    ]
fiber = monodromy_solve(fiber_cap_K, [a₀], [R[1,:]; r[1]], parameters = K)
fiber_points = solutions(fiber)


#
@polyvar f_var[1:N]
tracker = pathtracker([f - f_var; R * a[1:D-1] - r],
                        parameters=f_var,
                        start_parameters=f₀,
                        target_parameters=randn(N))
function move(p, tracker, fiber_points, f_SP)
    f₁ = f_SP(p)
    map(fiber_points) do p
        solution(track(tracker, p, details=:minimal, target_parameters = f₁))
    end
end
move_fiber(p) = move(p, tracker, fiber_points, f_SP)
# Monodromy by moving the linear space
@polyvar L[1:N+1]
X_cap_L = [
    transpose(L[1:N]) * f - L[N+1]; # moving linear condition
    Q[2:dimX, :] * f - q[2:dimX];   # fixed linear conditions
    R * a[1:D-1] - r                       # linear conditions for a
    ]
# monodromy
S = monodromy_solve(X_cap_L, fiber_points, [Q[1,:]; q[1]], parameters = L)
# extract solutions
sols = solutions(S)
f_sols = [f_SP(s) for s in sols]


X_cap_L_SP = StaticPolynomials.PolynomialSystem(X_cap_L)
[X_cap_L_SP([s; Q[1,:]; q[1]]) for in sols[1:2]]

#filter multiplicities
points = unique_points(f_sols)
sort(length.(multiplicities(f_sols)) )
degree_f = length(sols)/length(points)



test = map(s -> move_fiber(s), sols)
unique_points(vcat(test...))