restart;
with(LinearAlgebra):

# Define indeterminates

QParam := [q1,q2,q3,q4,q5,q6,q7,q8,q9,q10]:
PParam := Matrix([[p1],[p2],[p3],[p4],[p5],[p6],[p7],[p8],[p9],[p10],[p11],[p12],[p13],[p14],[p15]]):

# Special Fourier parameterization and inverse

F := Matrix([[1,1,1,1,1,1,1,1,1,1,1,1,1,1,1],
[1,-1/3,1,-1/3,1,-1/3,-1/3,1,-1/3,1,-1/3,1,-1/3,-1/3,1],
[1,5/9,1/9,1/9,11/27,7/27,1/9,1/9,-1/27,-5/27,-5/27,-1/27,-1/27,-5/27,-1/3],
[1,-1/3,-1/3,1/3,5/9,-1/3,-1/9,1/9,1/9,-1/3,1/3,1/9,-1/9,1/9,-1/3],
[1,-1/3,1,-1/3,1/9,-1/27,-1/27,1/9,-1/27,1/9,-1/27,-1/3,1/9,1/9,-1/3],
[1,3/7,-1/7,-1/7,1/21,-1/7,1/21,-1/3,-1/7,1/21,1/21,-1/21,-1/21,1/7,1/3],
[1,-1/3,-1/3,1/3,1/9,1/9,-1/9,-1/3,1/9,1/9,-1/9,-1/9,1/9,-1/9,1/3],
[1,1/3,-1/3,-1/3,-1/3,-1/21,-1/7,5/21,1/7,1/21,1/21,1/21,1/21,-1/21,-1/7],
[1,-1/3,-1/3,1/3,-1/3,1/21,1/7,5/21,-1/7,1/21,-1/21,1/21,-1/21,1/21,-1/7],
[1,-1/3,1,-1/3,-1/3,1/9,1/9,-1/3,1/9,-1/3,1/9,1/3,-1/9,-1/9,1/3]]):

FI := Matrix([
[1/256,3/256,27/256,9/128,27/256,21/128,9/64,21/128,21/128,9/128],
[3/128,-3/128,45/128,-9/64,-27/128,27/64,-9/32,21/64,-21/64,-9/64],
[3/256,9/256,9/256,-9/128,81/256,-9/128,-9/64,-21/128,-21/128,27/128],
[3/128,-3/128,9/128,9/64,-27/128,-9/64,9/32,-21/64,21/64,-9/64],
[9/256,27/256,99/256,45/128,27/256,9/128,9/64,-63/128,-63/128,-27/128],
[9/128,-9/128,63/128,-27/64,-9/128,-27/64,9/32,-9/64,9/64,9/64],
[9/64,-9/64,27/64,-9/32,-9/64,9/32,-9/16,-27/32,27/32,9/32],
[9/256,27/256,27/256,9/128,27/256,-63/128,-27/64,45/128,45/128,-27/128],
[9/64,-9/64,-9/64,9/32,-9/64,-27/32,9/16,27/32,-27/32,9/32],
[9/128,27/128,-45/128,-27/64,27/128,9/64,9/32,9/64,9/64,-27/64],
[9/128,-9/128,-45/128,27/64,-9/128,9/64,-9/32,9/64,-9/64,9/64],
[9/128,27/128,-9/128,9/64,-81/128,-9/64,-9/32,9/64,9/64,27/64],
[9/64,-9/64,-9/64,-9/32,27/64,-9/32,9/16,9/32,-9/32,-9/32],
[9/64,-9/64,-45/64,9/32,27/64,27/32,-9/16,-9/32,9/32,-9/32],
[3/128,9/128,-27/128,-9/64,-27/128,21/64,9/32,-9/64,-9/64,9/64]
]):

# List of polynomial parametrizations

P0 := [
d0^4*f0^5+3*d0^4*f1^5+9*d0^3*d1*f0^4*f1+3*d0^3*d1*f0^3*f1^2+3*d0^3*d1*f0^2*f1^3+9*d0^3*d1*f0*f1^4+24*d0^3*d1*f1^5+36*d0^2*d1^2*f0^3*f1^2+54*d0^2*d1^2*f0^2*f1^3+54*d0^2*d1^2*f0*f1^4+72*d0^2*d1^2*f1^5+9*d0*d1^3*f0^4*f1+39*d0*d1^3*f0^3*f1^2+99*d0*d1^3*f0^2*f1^3+189*d0*d1^3*f0*f1^4+96*d0*d1^3*f1^5+3*d1^4*f0^5+18*d1^4*f0^4*f1+42*d1^4*f0^3*f1^2+60*d1^4*f0^2*f1^3+72*d1^4*f0*f1^4+129*d1^4*f1^5,
6*d0^4*f0^4*f1+6*d0^4*f0*f1^4+12*d0^4*f1^5+6*d0^3*d1*f0^4*f1+66*d0^3*d1*f0^3*f1^2+18*d0^3*d1*f0^2*f1^3+102*d0^3*d1*f0*f1^4+96*d0^3*d1*f1^5+72*d0^2*d1^2*f0^3*f1^2+468*d0^2*d1^2*f0^2*f1^3+468*d0^2*d1^2*f0*f1^4+288*d0^2*d1^2*f1^5+6*d0*d1^3*f0^4*f1+138*d0*d1^3*f0^3*f1^2+738*d0*d1^3*f0^2*f1^3+1326*d0*d1^3*f0*f1^4+384*d0*d1^3*f1^5+30*d1^4*f0^4*f1+204*d1^4*f0^3*f1^2+504*d1^4*f0^2*f1^3+690*d1^4*f0*f1^4+516*d1^4*f1^5,
3*d0^4*f0^3*f1^2+3*d0^4*f0^2*f1^3+6*d0^4*f1^5+3*d0^3*d1*f0^5+15*d0^3*d1*f0^3*f1^2+51*d0^3*d1*f0^2*f1^3+18*d0^3*d1*f0*f1^4+57*d0^3*d1*f1^5+36*d0^2*d1^2*f0^4*f1+54*d0^2*d1^2*f0^3*f1^2+144*d0^2*d1^2*f0^2*f1^3+216*d0^2*d1^2*f0*f1^4+198*d0^2*d1^2*f1^5+3*d0*d1^3*f0^5+36*d0*d1^3*f0^4*f1+195*d0*d1^3*f0^3*f1^2+195*d0*d1^3*f0^2*f1^3+486*d0*d1^3*f0*f1^4+381*d0*d1^3*f1^5+6*d1^4*f0^5+36*d1^4*f0^4*f1+93*d1^4*f0^3*f1^2+255*d1^4*f0^2*f1^3+252*d1^4*f0*f1^4+330*d1^4*f1^5,
6*d0^4*f0^3*f1^2+12*d0^4*f0*f1^4+6*d0^4*f1^5+12*d0^3*d1*f0^4*f1+6*d0^3*d1*f0^3*f1^2+90*d0^3*d1*f0^2*f1^3+114*d0^3*d1*f0*f1^4+66*d0^3*d1*f1^5+144*d0^2*d1^2*f0^3*f1^2+288*d0^2*d1^2*f0^2*f1^3+612*d0^2*d1^2*f0*f1^4+252*d0^2*d1^2*f1^5+12*d0*d1^3*f0^4*f1+150*d0*d1^3*f0^3*f1^2+882*d0*d1^3*f0^2*f1^3+978*d0*d1^3*f0*f1^4+570*d0*d1^3*f1^5+24*d1^4*f0^4*f1+174*d1^4*f0^3*f1^2+468*d1^4*f0^2*f1^3+876*d1^4*f0*f1^4+402*d1^4*f1^5,
9*d0^4*f0^4*f1+9*d0^4*f0*f1^4+18*d0^4*f1^5+9*d0^3*d1*f0^5+18*d0^3*d1*f0^4*f1+63*d0^3*d1*f0^3*f1^2+63*d0^3*d1*f0^2*f1^3+108*d0^3*d1*f0*f1^4+171*d0^3*d1*f1^5+72*d0^2*d1^2*f0^4*f1+234*d0^2*d1^2*f0^3*f1^2+396*d0^2*d1^2*f0^2*f1^3+648*d0^2*d1^2*f0*f1^4+594*d0^2*d1^2*f1^5+9*d0*d1^3*f0^5+90*d0*d1^3*f0^4*f1+423*d0*d1^3*f0^3*f1^2+963*d0*d1^3*f0^2*f1^3+1260*d0*d1^3*f0*f1^4+1143*d0*d1^3*f1^5+18*d1^4*f0^5+135*d1^4*f0^4*f1+360*d1^4*f0^3*f1^2+522*d1^4*f0^2*f1^3+891*d1^4*f0*f1^4+990*d1^4*f1^5,
18*d0^4*f0^3*f1^2+18*d0^4*f0^2*f1^3+36*d0^4*f1^5+18*d0^3*d1*f0^4*f1+90*d0^3*d1*f0^3*f1^2+234*d0^3*d1*f0^2*f1^3+270*d0^3*d1*f0*f1^4+252*d0^3*d1*f1^5+36*d0^2*d1^2*f0^4*f1+288*d0^2*d1^2*f0^3*f1^2+1044*d0^2*d1^2*f0^2*f1^3+1764*d0^2*d1^2*f0*f1^4+756*d0^2*d1^2*f1^5+18*d0*d1^3*f0^4*f1+522*d0*d1^3*f0^3*f1^2+2322*d0*d1^3*f0^2*f1^3+3438*d0*d1^3*f0*f1^4+1476*d0*d1^3*f1^5+72*d1^4*f0^4*f1+522*d1^4*f0^3*f1^2+1566*d1^4*f0^2*f1^3+2304*d1^4*f0*f1^4+1368*d1^4*f1^5,
36*d0^4*f0^3*f1^2+72*d0^4*f0*f1^4+36*d0^4*f1^5+36*d0^3*d1*f0^4*f1+108*d0^3*d1*f0^3*f1^2+540*d0^3*d1*f0^2*f1^3+612*d0^3*d1*f0*f1^4+432*d0^3*d1*f1^5+36*d0^2*d1^2*f0^4*f1+648*d0^2*d1^2*f0^3*f1^2+2160*d0^2*d1^2*f0^2*f1^3+3312*d0^2*d1^2*f0*f1^4+1620*d0^2*d1^2*f1^5+108*d0*d1^3*f0^4*f1+1116*d0*d1^3*f0^3*f1^2+4428*d0*d1^3*f0^2*f1^3+6804*d0*d1^3*f0*f1^4+3096*d0*d1^3*f1^5+108*d1^4*f0^4*f1+972*d1^4*f0^3*f1^2+3240*d1^4*f0^2*f1^3+4752*d1^4*f0*f1^4+2592*d1^4*f1^5,
9*d0^4*f0^3*f1^2+9*d0^4*f0^2*f1^3+18*d0^4*f1^5+27*d0^3*d1*f0^4*f1+45*d0^3*d1*f0^3*f1^2+81*d0^3*d1*f0^2*f1^3+135*d0^3*d1*f0*f1^4+144*d0^3*d1*f1^5+18*d0^2*d1^2*f0^5+54*d0^2*d1^2*f0^4*f1+216*d0^2*d1^2*f0^3*f1^2+414*d0^2*d1^2*f0^2*f1^3+648*d0^2*d1^2*f0*f1^4+594*d0^2*d1^2*f1^5+135*d0*d1^3*f0^4*f1+477*d0*d1^3*f0^3*f1^2+837*d0*d1^3*f0^2*f1^3+1215*d0*d1^3*f0*f1^4+1224*d0*d1^3*f1^5+18*d1^4*f0^5+108*d1^4*f0^4*f1+333*d1^4*f0^3*f1^2+603*d1^4*f0^2*f1^3+918*d1^4*f0*f1^4+936*d1^4*f1^5,
72*d0^4*f0^2*f1^3+36*d0^4*f0*f1^4+36*d0^4*f1^5+180*d0^3*d1*f0^3*f1^2+396*d0^3*d1*f0^2*f1^3+828*d0^3*d1*f0*f1^4+324*d0^3*d1*f1^5+108*d0^2*d1^2*f0^4*f1+504*d0^2*d1^2*f0^3*f1^2+2232*d0^2*d1^2*f0^2*f1^3+3312*d0^2*d1^2*f0*f1^4+1620*d0^2*d1^2*f1^5+72*d0*d1^3*f0^4*f1+1332*d0*d1^3*f0^3*f1^2+4428*d0*d1^3*f0^2*f1^3+6300*d0*d1^3*f0*f1^4+3420*d0*d1^3*f1^5+108*d1^4*f0^4*f1+864*d1^4*f0^3*f1^2+3240*d1^4*f0^2*f1^3+5076*d1^4*f0*f1^4+2376*d1^4*f1^5,
36*d0^4*f0^2*f1^3+18*d0^4*f0*f1^4+18*d0^4*f1^5+18*d0^3*d1*f0^4*f1+90*d0^3*d1*f0^3*f1^2+234*d0^3*d1*f0^2*f1^3+270*d0^3*d1*f0*f1^4+252*d0^3*d1*f1^5+18*d0^2*d1^2*f0^5+126*d0^2*d1^2*f0^4*f1+504*d0^2*d1^2*f0^3*f1^2+756*d0^2*d1^2*f0^2*f1^3+1242*d0^2*d1^2*f0*f1^4+1242*d0^2*d1^2*f1^5+36*d0*d1^3*f0^5+342*d0*d1^3*f0^4*f1+882*d0*d1^3*f0^3*f1^2+1530*d0*d1^3*f0^2*f1^3+2466*d0*d1^3*f0*f1^4+2520*d0*d1^3*f1^5+18*d1^4*f0^5+162*d1^4*f0^4*f1+684*d1^4*f0^3*f1^2+1332*d1^4*f0^2*f1^3+1836*d1^4*f0*f1^4+1800*d1^4*f1^5,
18*d0^4*f0^2*f1^3+54*d0^4*f0*f1^4+54*d0^3*d1*f0^3*f1^2+234*d0^3*d1*f0^2*f1^3+450*d0^3*d1*f0*f1^4+126*d0^3*d1*f1^5+36*d0^2*d1^2*f0^4*f1+288*d0^2*d1^2*f0^3*f1^2+1152*d0^2*d1^2*f0^2*f1^3+1548*d0^2*d1^2*f0*f1^4+864*d0^2*d1^2*f1^5+72*d0*d1^3*f0^4*f1+702*d0*d1^3*f0^3*f1^2+2106*d0*d1^3*f0^2*f1^3+3114*d0*d1^3*f0*f1^4+1782*d0*d1^3*f1^5+36*d1^4*f0^4*f1+396*d1^4*f0^3*f1^2+1674*d1^4*f0^2*f1^3+2610*d1^4*f0*f1^4+1116*d1^4*f1^5,
18*d0^4*f0^3*f1^2+36*d0^4*f0*f1^4+18*d0^4*f1^5+36*d0^3*d1*f0^4*f1+108*d0^3*d1*f0^3*f1^2+144*d0^3*d1*f0^2*f1^3+324*d0^3*d1*f0*f1^4+252*d0^3*d1*f1^5+18*d0^2*d1^2*f0^5+144*d0^2*d1^2*f0^4*f1+414*d0^2*d1^2*f0^3*f1^2+882*d0^2*d1^2*f0^2*f1^3+1188*d0^2*d1^2*f0*f1^4+1242*d0^2*d1^2*f1^5+36*d0*d1^3*f0^5+252*d0*d1^3*f0^4*f1+936*d0*d1^3*f0^3*f1^2+1692*d0*d1^3*f0^2*f1^3+2340*d0*d1^3*f0*f1^4+2520*d0*d1^3*f1^5+18*d1^4*f0^5+216*d1^4*f0^4*f1+684*d1^4*f0^3*f1^2+1170*d1^4*f0^2*f1^3+1944*d1^4*f0*f1^4+1800*d1^4*f1^5,
72*d0^4*f0^2*f1^3+36*d0^4*f0*f1^4+36*d0^4*f1^5+144*d0^3*d1*f0^3*f1^2+504*d0^3*d1*f0^2*f1^3+720*d0^3*d1*f0*f1^4+360*d0^3*d1*f1^5+72*d0^2*d1^2*f0^4*f1+612*d0^2*d1^2*f0^3*f1^2+2124*d0^2*d1^2*f0^2*f1^3+3348*d0^2*d1^2*f0*f1^4+1620*d0^2*d1^2*f1^5+144*d0*d1^3*f0^4*f1+1224*d0*d1^3*f0^3*f1^2+4320*d0*d1^3*f0^2*f1^3+6552*d0*d1^3*f0*f1^4+3312*d0*d1^3*f1^5+72*d1^4*f0^4*f1+900*d1^4*f0^3*f1^2+3348*d1^4*f0^2*f1^3+4896*d1^4*f0*f1^4+2448*d1^4*f1^5,
36*d0^4*f0^2*f1^3+108*d0^4*f0*f1^4+72*d0^3*d1*f0^3*f1^2+576*d0^3*d1*f0^2*f1^3+792*d0^3*d1*f0*f1^4+288*d0^3*d1*f1^5+36*d0^2*d1^2*f0^4*f1+684*d0^2*d1^2*f0^3*f1^2+2196*d0^2*d1^2*f0^2*f1^3+3132*d0^2*d1^2*f0*f1^4+1728*d0^2*d1^2*f1^5+216*d0*d1^3*f0^4*f1+1296*d0*d1^3*f0^3*f1^2+4104*d0*d1^3*f0^2*f1^3+6480*d0*d1^3*f0*f1^4+3456*d0*d1^3*f1^5+36*d1^4*f0^4*f1+828*d1^4*f0^3*f1^2+3456*d1^4*f0^2*f1^3+5040*d1^4*f0*f1^4+2304*d1^4*f1^5,
6*d0^4*f0^2*f1^3+18*d0^4*f0*f1^4+36*d0^3*d1*f0^3*f1^2+72*d0^3*d1*f0^2*f1^3+108*d0^3*d1*f0*f1^4+72*d0^3*d1*f1^5+54*d0^2*d1^2*f0^4*f1+162*d0^2*d1^2*f0^3*f1^2+270*d0^2*d1^2*f0^2*f1^3+378*d0^2*d1^2*f0*f1^4+432*d0^2*d1^2*f1^5+24*d0*d1^3*f0^5+108*d0*d1^3*f0^4*f1+288*d0*d1^3*f0^3*f1^2+516*d0*d1^3*f0^2*f1^3+792*d0*d1^3*f0*f1^4+864*d0*d1^3*f1^5+54*d1^4*f0^4*f1+234*d1^4*f0^3*f1^2+432*d1^4*f0^2*f1^3+648*d1^4*f0*f1^4+576*d1^4*f1^5]:

# Substitutions based on the model

P := P0:
P := subs(d0 = 1-3*d1, P):
P := subs(f0 = 1-3*f1, 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([
q7*q9-q8*q10,
q5*q9-q7*q10,
q4*q9-q6*q10,
q7*q8-q6*q9,
q5*q8-q6*q10,
q4*q8-q3*q9,
q7^2-q6*q10,
q6*q7-q3*q9,
q5*q7-q4*q10,
q4*q7-q3*q10,
q6^2-q3*q8,
q5*q6-q3*q10,
q4*q6-q3*q7,
q2*q6-q1*q10,
q4*q5-q2*q9,
q3*q5-q2*q8,
q4^2-q2*q8,
q3*q4-q1*q9,
q3^2-q1*q8,
q2*q3-q1*q5,
q5^3-q2*q10^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: