# Code for Example 3.4 in "Gibbs Manifolds" by Dmitrii Pavlov, Bernd Sturmfels and Simon Telen
# Julia v1.8.3
# November 22, 2022
using Oscar
# Make a 4x4 Hankel matrix from a length 7 vector
function hankel(q)
[q[1:4]';q[2:5]';q[3:6]';q[4:7]']
end
# Generate enough samples on the Gibbs manifold to interpolate in degree d
d = 6
nsamples = binomial(9+d,d)+200
sample_matrices = [exp(hankel([0;2*pi*im*rand(6)])) for i = 1:nsamples]
sample_vectors = [[M[i,j] for i = 1:4 for j = i:4] for M in sample_matrices]
# Define a polynomial ring for the equations of the Gibbs variety
x_strings = ["x$i$j" for i = 1:4 for j = i:4]
R, x = PolynomialRing(QQ,x_strings)
# Make a Vandermonde matrix by evaluating all monomials of degree d at all sample vectors
exp_vecs = collect(exponent_vectors(sum(x)^d))
vdm_matrix = transpose(hcat([[prod(sv.^ev) for ev in exp_vecs] for sv in sample_vectors]...))
# We normalize the rows for numerical stability, and compute the kernel
for i = 1:size(vdm_matrix,1)
vdm_matrix[i,:] = vdm_matrix[i,:]/norm(vdm_matrix[i,:])
end
N = nullspace(vdm_matrix)
# We rationalize the entries of a normalized basis vector of the kernel
N = N/N[findfirst(k->abs(k)>1e-8, N)]
N[findall(abs.(N) .<1e-10)] .= 0
N = real(N)
N = rationalize.(N,tol=1e-7)
# We check that the rationalized vector is still in the kernel
norm(vdm_matrix*N)
# We construct the interpolating polynomial
mons = [prod(x.^ev) for ev in exp_vecs]
pol = (transpose(N)*mons)[1]
length(terms(pol))
# We compute the Newton polytope and its f-vector
P = newton_polytope(pol)
f_vector(P)