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/9,1/9,-1/3,1,1/9,-1/3,1,1/9,1/9,-1/3,1/9,-1/3,1,1/9],
[1,-1/3,-1/3,1/3,1,-1/3,1/3,1,-1/3,-1/3,1/3,-1/3,1/3,1,-1/3],
[1,13/21,5/21,5/21,-1/7,-1/7,-1/7,5/21,5/21,1/21,1/21,-1/7,-1/7,-1/3,-1/3],
[1,1/9,-1/3,-1/9,-1/3,1/9,1/3,1/3,-1/9,1/9,-1/9,-1/9,1/9,-1/3,1/9],
[1,-11/69,13/69,-7/69,13/69,5/69,-5/23,5/69,-1/23,-7/69,5/69,5/69,1/69,-1/3,1/69],
[1,-1/3,-1/3,1/3,1,-1/3,1/3,-1/3,1/9,1/9,-1/9,1/9,-1/9,-1/3,1/9],
[1,5/9,1/9,1/9,-1/3,-1/3,-1/3,-1/3,-1/3,-1/9,-1/9,1/9,1/9,1/3,1/3],
[1,1/9,-1/3,-1/9,-1/3,1/9,1/3,-1/3,1/9,-1/9,1/9,1/9,-1/9,1/3,-1/9],
[1,-1/3,5/21,1/21,-1/3,1/21,-1/7,-1/3,1/21,1/7,-1/21,-1/7,1/21,1/3,-1/21]]):

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

# List of polynomial parametrizations

P0 := [
c0^3*f0^5+3*c0^3*f1^5+6*c0^2*c1*f0^4*f1+3*c0^2*c1*f0^3*f1^2+3*c0^2*c1*f0^2*f1^3+6*c0^2*c1*f0*f1^4+18*c0^2*c1*f1^5+6*c0*c1^2*f0^4*f1+21*c0*c1^2*f0^3*f1^2+3*c0*c1^2*f0^2*f1^3+42*c0*c1^2*f0*f1^4+36*c0*c1^2*f1^5+3*c1^3*f0^5+12*c1^3*f0^4*f1+12*c1^3*f0^3*f1^2+6*c1^3*f0^2*f1^3+24*c1^3*f0*f1^4+51*c1^3*f1^5,
9*c0^3*f0^4*f1+9*c0^3*f0*f1^4+18*c0^3*f1^5+63*c0^2*c1*f0^3*f1^2+63*c0^2*c1*f0^2*f1^3+90*c0^2*c1*f0*f1^4+108*c0^2*c1*f1^5+63*c0*c1^2*f0^3*f1^2+333*c0*c1^2*f0^2*f1^3+360*c0*c1^2*f0*f1^4+216*c0*c1^2*f1^5+27*c1^3*f0^4*f1+126*c1^3*f0^3*f1^2+216*c1^3*f0^2*f1^3+297*c1^3*f0*f1^4+306*c1^3*f1^5,
9*c0^3*f0^3*f1^2+9*c0^3*f0^2*f1^3+18*c0^3*f1^5+9*c0^2*c1*f0^4*f1+18*c0^2*c1*f0^3*f1^2+126*c0^2*c1*f0^2*f1^3+63*c0^2*c1*f0*f1^4+108*c0^2*c1*f1^5+9*c0*c1^2*f0^4*f1+126*c0*c1^2*f0^3*f1^2+180*c0*c1^2*f0^2*f1^3+441*c0*c1^2*f0*f1^4+216*c0*c1^2*f1^5+18*c1^3*f0^4*f1+99*c1^3*f0^3*f1^2+297*c1^3*f0^2*f1^3+252*c1^3*f0*f1^4+306*c1^3*f1^5,
18*c0^3*f0^3*f1^2+36*c0^3*f0*f1^4+18*c0^3*f1^5+36*c0^2*c1*f0^3*f1^2+198*c0^2*c1*f0^2*f1^3+306*c0^2*c1*f0*f1^4+108*c0^2*c1*f1^5+36*c0*c1^2*f0^3*f1^2+630*c0*c1^2*f0^2*f1^3+1062*c0*c1^2*f0*f1^4+216*c0*c1^2*f1^5+126*c1^3*f0^3*f1^2+540*c1^3*f0^2*f1^3+972*c1^3*f0*f1^4+306*c1^3*f1^5,
3*c0^3*f0^3*f1^2+3*c0^3*f0^2*f1^3+6*c0^3*f1^5+3*c0^2*c1*f0^5+6*c0^2*c1*f0^4*f1+18*c0^2*c1*f0^3*f1^2+6*c0^2*c1*f0^2*f1^3+30*c0^2*c1*f0*f1^4+45*c0^2*c1*f1^5+3*c0*c1^2*f0^5+42*c0*c1^2*f0^4*f1+36*c0*c1^2*f0^3*f1^2+6*c0*c1^2*f0^2*f1^3+102*c0*c1^2*f0*f1^4+135*c0*c1^2*f1^5+6*c1^3*f0^5+24*c1^3*f0^4*f1+51*c1^3*f0^3*f1^2+21*c1^3*f0^2*f1^3+84*c1^3*f0*f1^4+138*c1^3*f1^5,
36*c0^3*f0^2*f1^3+18*c0^3*f0*f1^4+18*c0^3*f1^5+18*c0^2*c1*f0^4*f1+54*c0^2*c1*f0^3*f1^2+162*c0^2*c1*f0^2*f1^3+252*c0^2*c1*f0*f1^4+162*c0^2*c1*f1^5+18*c0*c1^2*f0^4*f1+270*c0*c1^2*f0^3*f1^2+486*c0*c1^2*f0^2*f1^3+576*c0*c1^2*f0*f1^4+594*c0*c1^2*f1^5+36*c1^3*f0^4*f1+180*c1^3*f0^3*f1^2+540*c1^3*f0^2*f1^3+666*c1^3*f0*f1^4+522*c1^3*f1^5,
6*c0^3*f0^2*f1^3+18*c0^3*f0*f1^4+18*c0^2*c1*f0^3*f1^2+36*c0^2*c1*f0^2*f1^3+144*c0^2*c1*f0*f1^4+18*c0^2*c1*f1^5+18*c0*c1^2*f0^3*f1^2+252*c0*c1^2*f0^2*f1^3+252*c0*c1^2*f0*f1^4+126*c0*c1^2*f1^5+36*c1^3*f0^3*f1^2+162*c1^3*f0^2*f1^3+378*c1^3*f0*f1^4+72*c1^3*f1^5,
6*c0^3*f0^4*f1+6*c0^3*f0*f1^4+12*c0^3*f1^5+6*c0^2*c1*f0^5+18*c0^2*c1*f0^4*f1+30*c0^2*c1*f0^3*f1^2+6*c0^2*c1*f0^2*f1^3+66*c0^2*c1*f0*f1^4+90*c0^2*c1*f1^5+6*c0*c1^2*f0^5+54*c0*c1^2*f0^4*f1+102*c0*c1^2*f0^3*f1^2+42*c0*c1^2*f0^2*f1^3+174*c0*c1^2*f0*f1^4+270*c0*c1^2*f1^5+12*c1^3*f0^5+66*c1^3*f0^4*f1+84*c1^3*f0^3*f1^2+24*c1^3*f0^2*f1^3+186*c1^3*f0*f1^4+276*c1^3*f1^5,
18*c0^3*f0^3*f1^2+18*c0^3*f0^2*f1^3+36*c0^3*f1^5+18*c0^2*c1*f0^4*f1+72*c0^2*c1*f0^3*f1^2+144*c0^2*c1*f0^2*f1^3+234*c0^2*c1*f0*f1^4+180*c0^2*c1*f1^5+18*c0*c1^2*f0^4*f1+180*c0*c1^2*f0^3*f1^2+576*c0*c1^2*f0^2*f1^3+666*c0*c1^2*f0*f1^4+504*c0*c1^2*f1^5+36*c1^3*f0^4*f1+234*c1^3*f0^3*f1^2+486*c1^3*f0^2*f1^3+612*c1^3*f0*f1^4+576*c1^3*f1^5,
36*c0^3*f0^3*f1^2+72*c0^3*f0*f1^4+36*c0^3*f1^5+36*c0^2*c1*f0^4*f1+108*c0^2*c1*f0^3*f1^2+360*c0^2*c1*f0^2*f1^3+432*c0^2*c1*f0*f1^4+360*c0^2*c1*f1^5+36*c0*c1^2*f0^4*f1+432*c0*c1^2*f0^3*f1^2+1116*c0*c1^2*f0^2*f1^3+1188*c0*c1^2*f0*f1^4+1116*c0*c1^2*f1^5+72*c1^3*f0^4*f1+432*c1^3*f0^3*f1^2+972*c1^3*f0^2*f1^3+1332*c1^3*f0*f1^4+1080*c1^3*f1^5,
72*c0^3*f0^2*f1^3+36*c0^3*f0*f1^4+36*c0^3*f1^5+72*c0^2*c1*f0^3*f1^2+360*c0^2*c1*f0^2*f1^3+684*c0^2*c1*f0*f1^4+180*c0^2*c1*f1^5+180*c0*c1^2*f0^3*f1^2+1116*c0*c1^2*f0^2*f1^3+2088*c0*c1^2*f0*f1^4+504*c0*c1^2*f1^5+180*c1^3*f0^3*f1^2+1188*c1^3*f0^2*f1^3+1944*c1^3*f0*f1^4+576*c1^3*f1^5,
72*c0^3*f0^2*f1^3+36*c0^3*f0*f1^4+36*c0^3*f1^5+180*c0^2*c1*f0^3*f1^2+324*c0^2*c1*f0^2*f1^3+432*c0^2*c1*f0*f1^4+360*c0^2*c1*f1^5+108*c0*c1^2*f0^4*f1+396*c0*c1^2*f0^3*f1^2+972*c0*c1^2*f0^2*f1^3+1296*c0*c1^2*f0*f1^4+1116*c0*c1^2*f1^5+36*c1^3*f0^4*f1+432*c1^3*f0^3*f1^2+1080*c1^3*f0^2*f1^3+1260*c1^3*f0*f1^4+1080*c1^3*f1^5,
36*c0^3*f0^2*f1^3+108*c0^3*f0*f1^4+36*c0^2*c1*f0^3*f1^2+432*c0^2*c1*f0^2*f1^3+648*c0^2*c1*f0*f1^4+180*c0^2*c1*f1^5+252*c0*c1^2*f0^3*f1^2+1080*c0*c1^2*f0^2*f1^3+1944*c0*c1^2*f0*f1^4+612*c0*c1^2*f1^5+144*c1^3*f0^3*f1^2+1188*c1^3*f0^2*f1^3+2052*c1^3*f0*f1^4+504*c1^3*f1^5,
6*c0^3*f0^3*f1^2+12*c0^3*f0*f1^4+6*c0^3*f1^5+24*c0^2*c1*f0^4*f1+30*c0^2*c1*f0^3*f1^2+12*c0^2*c1*f0^2*f1^3+60*c0^2*c1*f0*f1^4+90*c0^2*c1*f1^5+18*c0*c1^2*f0^5+60*c0*c1^2*f0^4*f1+84*c0*c1^2*f0^3*f1^2+30*c0*c1^2*f0^2*f1^3+168*c0*c1^2*f0*f1^4+288*c0*c1^2*f1^5+6*c1^3*f0^5+60*c1^3*f0^4*f1+96*c1^3*f0^3*f1^2+30*c1^3*f0^2*f1^3+192*c1^3*f0*f1^4+264*c1^3*f1^5,
18*c0^3*f0^2*f1^3+54*c0^3*f0*f1^4+72*c0^2*c1*f0^3*f1^2+198*c0^2*c1*f0^2*f1^3+198*c0^2*c1*f0*f1^4+180*c0^2*c1*f1^5+54*c0*c1^2*f0^4*f1+234*c0*c1^2*f0^3*f1^2+468*c0*c1^2*f0^2*f1^3+576*c0*c1^2*f0*f1^4+612*c0*c1^2*f1^5+18*c1^3*f0^4*f1+198*c1^3*f0^3*f1^2+540*c1^3*f0^2*f1^3+684*c1^3*f0*f1^4+504*c1^3*f1^5]:

# Substitutions based on the model

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