restart;
with(LinearAlgebra):

# Define indeterminates

QParam := [q1,q2,q3,q4,q5,q6,q7]:
PParam := Matrix([[p1],[p2],[p3],[p4],[p5],[p6],[p7],[p8],[p9]]):

# Special Fourier parameterization and inverse

F := Matrix([
[1,1,1,1,1,1,1,1,1],
[1,-1/3,1,-1/3,1,-1/3,-1/3,1,-1/3],
[1,7/15,-1/15,-1/15,1/5,1/5,-1/15,-1/3,-1/3],
[1,-1/3,-1/3,1/3,1/3,-1/3,0,-1/3,1/3],
[1,-1/3,1,-1/3,-1/3,1/9,1/9,-1/3,1/9],
[1,1/3,-1/3,-1/3,-1/3,-1/3,0,1/3,1/3],
[1,-1/3,-1/3,1/3,-1/3,1/3,0,1/3,-1/3]]):

FI := Matrix([
[1/64,3/64,15/64,3/16,9/64,3/16,3/16],
[3/32,-3/32,21/32,-3/8,-9/32,3/8,-3/8],
[3/64,9/64,-3/64,-3/16,27/64,-3/16,-3/16],
[3/32,-3/32,-3/32,3/8,-9/32,-3/8,3/8],
[3/32,9/32,9/32,3/8,-9/32,-3/8,-3/8],
[3/32,-3/32,9/32,-3/8,3/32,-3/8,3/8],
[3/8,-3/8,-3/8,0,3/8,0,0],
[3/32,9/32,-15/32,-3/8,-9/32,3/8,3/8],
[3/32,-3/32,-15/32,3/8,3/32,3/8,-3/8]]):

# List of polynomial parametrizations

P0 := [
c0^3*e0^4+3*c0^3*e1^4+6*c0^2*c1*e0^3*e1+6*c0^2*c1*e0^2*e1^2+6*c0^2*c1*e0*e1^3+18*c0^2*c1*e1^4+6*c0*c1^2*e0^3*e1+24*c0*c1^2*e0^2*e1^2+42*c0*c1^2*e0*e1^3+36*c0*c1^2*e1^4+3*c1^3*e0^4+12*c1^3*e0^3*e1+18*c1^3*e0^2*e1^2+24*c1^3*e0*e1^3+51*c1^3*e1^4,
6*c0^3*e0^3*e1+6*c0^3*e0*e1^3+12*c0^3*e1^4+6*c0^2*c1*e0^3*e1+60*c0^2*c1*e0^2*e1^2+78*c0^2*c1*e0*e1^3+72*c0^2*c1*e1^4+6*c0*c1^2*e0^3*e1+132*c0*c1^2*e0^2*e1^2+366*c0*c1^2*e0*e1^3+144*c0*c1^2*e1^4+30*c1^3*e0^3*e1+144*c1^3*e0^2*e1^2+270*c1^3*e0*e1^3+204*c1^3*e1^4,
6*c0^3*e0^2*e1^2+6*c0^3*e1^4+3*c0^2*c1*e0^4+6*c0^2*c1*e0^3*e1+24*c0^2*c1*e0^2*e1^2+30*c0^2*c1*e0*e1^3+45*c0^2*c1*e1^4+3*c0*c1^2*e0^4+42*c0*c1^2*e0^3*e1+42*c0*c1^2*e0^2*e1^2+102*c0*c1^2*e0*e1^3+135*c0*c1^2*e1^4+6*c1^3*e0^4+24*c1^3*e0^3*e1+72*c1^3*e0^2*e1^2+84*c1^3*e0*e1^3+138*c1^3*e1^4,
6*c0^3*e0^2*e1^2+12*c0^3*e0*e1^3+6*c0^3*e1^4+12*c0^2*c1*e0^3*e1+30*c0^2*c1*e0^2*e1^2+120*c0^2*c1*e0*e1^3+54*c0^2*c1*e1^4+12*c0*c1^2*e0^3*e1+174*c0*c1^2*e0^2*e1^2+264*c0*c1^2*e0*e1^3+198*c0*c1^2*e1^4+24*c1^3*e0^3*e1+126*c1^3*e0^2*e1^2+324*c1^3*e0*e1^3+174*c1^3*e1^4,
6*c0^3*e0^3*e1+6*c0^3*e0*e1^3+12*c0^3*e1^4+6*c0^2*c1*e0^4+18*c0^2*c1*e0^3*e1+36*c0^2*c1*e0^2*e1^2+66*c0^2*c1*e0*e1^3+90*c0^2*c1*e1^4+6*c0*c1^2*e0^4+54*c0*c1^2*e0^3*e1+144*c0*c1^2*e0^2*e1^2+174*c0*c1^2*e0*e1^3+270*c0*c1^2*e1^4+12*c1^3*e0^4+66*c1^3*e0^3*e1+108*c1^3*e0^2*e1^2+186*c1^3*e0*e1^3+276*c1^3*e1^4,
12*c0^3*e0^2*e1^2+12*c0^3*e1^4+12*c0^2*c1*e0^3*e1+36*c0^2*c1*e0^2*e1^2+108*c0^2*c1*e0*e1^3+60*c0^2*c1*e1^4+12*c0*c1^2*e0^3*e1+144*c0*c1^2*e0^2*e1^2+324*c0*c1^2*e0*e1^3+168*c0*c1^2*e1^4+24*c1^3*e0^3*e1+144*c1^3*e0^2*e1^2+288*c1^3*e0*e1^3+192*c1^3*e1^4,
24*c0^3*e0^2*e1^2+48*c0^3*e0*e1^3+24*c0^3*e1^4+24*c0^2*c1*e0^3*e1+192*c0^2*c1*e0^2*e1^2+408*c0^2*c1*e0*e1^3+240*c0^2*c1*e1^4+96*c0*c1^2*e0^3*e1+552*c0*c1^2*e0^2*e1^2+1200*c0*c1^2*e0*e1^3+744*c0*c1^2*e1^4+72*c1^3*e0^3*e1+576*c1^3*e0^2*e1^2+1224*c1^3*e0*e1^3+720*c1^3*e1^4,
6*c0^3*e0^2*e1^2+12*c0^3*e0*e1^3+6*c0^3*e1^4+24*c0^2*c1*e0^3*e1+42*c0^2*c1*e0^2*e1^2+60*c0^2*c1*e0*e1^3+90*c0^2*c1*e1^4+18*c0*c1^2*e0^4+60*c0*c1^2*e0^3*e1+114*c0*c1^2*e0^2*e1^2+168*c0*c1^2*e0*e1^3+288*c0*c1^2*e1^4+6*c1^3*e0^4+60*c1^3*e0^3*e1+126*c1^3*e0^2*e1^2+192*c1^3*e0*e1^3+264*c1^3*e1^4,
24*c0^3*e0*e1^3+60*c0^2*c1*e0^2*e1^2+96*c0^2*c1*e0*e1^3+60*c0^2*c1*e1^4+36*c0*c1^2*e0^3*e1+132*c0*c1^2*e0^2*e1^2+276*c0*c1^2*e0*e1^3+204*c0*c1^2*e1^4+12*c1^3*e0^3*e1+144*c1^3*e0^2*e1^2+324*c1^3*e0*e1^3+168*c1^3*e1^4
]:

# Substitutions based on the model

P := P0:
P := subs(c0 = 1-3*c1, P):
P := subs(e0 = 1-3*e1, P):

# Check that the polynomial parametrization lies in the probability simplex

suma := 0:
for i from 1 to nops(P) do
suma := suma + P[i]:
od:
normal(expand(suma));

# Ideal of Invariants in Fourier coordinates

Invariants := Matrix([
q5*q6-q4*q7,
q4*q6-q3*q7,
q4^2-q3*q5,
q2*q3-q1*q5,
q4^2*q5-q2*q7^2,
q4^3-q2*q6*q7,
q3*q4^2-q2*q6^2,
q3*q4^2-q1*q7^2,
q3^2*q4-q1*q6*q7,
q3^3-q1*q6^2]):

# Ideal of Invariants in probability coordinates

Fourier := MatrixMatrixMultiply(F,PParam):

PInvariants := Invariants:
for i from 1 to nops(QParam) do
PInvariants := subs(QParam[i] = Fourier[i, 1], PInvariants):
od:

# Evaluation of Invariants at the polynomial/rational parametrization

num := op(PInvariants[1,1..-1])[1]:
for j from 1 to num do
coordpoly := PInvariants[1, j]:
for i from 1 to op(PParam[1..-1,1])[1] do
coordpoly := subs(PParam[i, 1] = P0[i], coordpoly):
od:
coordpoly :=expand(coordpoly):
lprint(j,coordpoly);
od: