//Author: Daniel Hughes

function MaxPS(g)
PN:=parabolicNames[g];
n:=#PN;
pn:=PN;
for i:=1 to n-1 do
Pi:=paraData(g,PN[i]);
for j:=i+1 to n do
Pj:=paraData(g,PN[j]);
if Pi meet Pj eq Pi then Exclude(~pn,PN[i]);
elif Pi meet Pj eq Pj then Exclude(~pn,PN[j]);
end if;
end for;
end for;
return pn;
end function;

function OrdMaxPS(g)
MPS:=MaxPS(g);Orders:=[];
for i:=1 to #MPS do
Pi:=paraData(g,MPS[i]);
Orders[i]:=Order(Pi);
end for;
return Orders;
end function;

function QOrdMaxPS(g)
MPS:=MaxPS(g);Orders:=[];G:=ShephardTodd(g);
for i:=1 to #MPS do
Pi:=paraData(g,MPS[i]);
Orders[i]:=Order(G)/Order(Pi);
end for;
return Orders;
end function;

HermTranspose := function(B)
T:=Transpose(B);
T:=Parent(T)![ComplexConjugate(z):z in Eltseq(T)];
return T;
end function;

ip := function(u,v)
if Type(BaseRing(u)) in {FldCyc,FldCom,FldRat,FldRe}
then return InnerProduct(u,v);
else return "The inner product is not defined over this field!";
end if;
end function;

function inp(a,b)
innerp:=&+[a[i]*ComplexConjugate(b[i]):i in [1..Ncols(a)]];
return(innerp);
end function;

CanonicalGramian := function(V)
F:=BaseRing(V); n:=Ncols(V);
dep:=NullspaceOfTranspose(V); B:=Basis(dep);
L:=[];
for j in [1..#B] do
v:=B[j]; prj:=0*v;
for k in [1..j-1] do
prj:=prj+inp(v,L[k])/inp(L[k],L[k])*L[k];
end for;
Append(~L,v-prj);
end for;
P:=Matrix(F,n,n,[]);
for j in [1..#L] do
v:=L[j]; InnerP:=inp(v,v); v:=Matrix(F,n,1,[v[j]:j in [1..n]]);
P:=P+(1/InnerP)*v*HermTranspose(v);
end for;
I:=Identity(Parent(P)); // I:=Matrix(F,n,n,[<j,j,1>: j in [1..n]]);
P:=I-P;
return(P);
end function;


function tCalc(g)
MPS:=MaxPS(g);tvals:=[];Avals:=[* *];
for i:=1 to #MPS do

P,Q,H,N,roots,rho,J,Xi,G:=paraData(g,MPS[i]);
a:=&meet[Eigenspace(gen,1):gen in P].1;
T:=RightTransversal(G,P);
v:=[a^g:g in T];
V:=Transpose(Matrix(v));
A:=CanonicalGramian(V);
d:=Nrows(V);

dif:=0; t:=0;
while dif eq 0 do
t:=t+1;
if BaseRing(A) eq Rationals() then
c:=&*[(1+i)/(d+i): i in [0..2*(t-1) by 2]];
else
c:=Factorial(d-1)*Factorial(t)/Factorial(d+t-1);
end if;

n:=Ncols(A); sum:=0;
for i:=1 to n do
for j:=1 to n do
sum:=sum + (A[i,j]*ComplexConjugate(A[i,j]))^(t);
end for;
end for;

D:=Diagonal(A); diagsum:=0;
for i:=1 to n do
diagsum:=diagsum + D[i]^t;
end for;
diagsum:=diagsum^2;

dif:=sum-c*diagsum;

end while;

Append(~tvals,t-1);Append(~Avals,A);
end for;

return tvals, Avals;
end function;

function tCalci(g,i)
MPS:=MaxPS(g);tvals:=[];Avals:=[];

P,Q,H,N,roots,rho,J,Xi,G:=paraData(g,MPS[i]);
a:=&meet[Eigenspace(gen,1):gen in P].1;
T:=RightTransversal(G,P);
v:=[a^g:g in T];
V:=Transpose(Matrix(v));
A:=CanonicalGramian(V);
d:=Nrows(V);

dif:=0; t:=0;
while dif eq 0 do
t:=t+1;
if BaseRing(A) eq Rationals() then
c:=&*[(1+i)/(d+i): i in [0..2*(t-1) by 2]];
else c:=Factorial(d-1)*Factorial(t)/Factorial(d+t-1);

n:=Ncols(A); sum:=0;
for i:=1 to n do
for j:=1 to n do
sum:=sum + (A[i,j]*ComplexConjugate(A[i,j]))^(t);
end for;
end for;

D:=Diagonal(A); diagsum:=0;
for i:=1 to n do
diagsum:=diagsum + D[i]^t;
end for;
diagsum:=diagsum^2;

dif:=sum-c*diagsum;

end while;

return t-1, A;
end function;

load "ComputeSymmetry.magma";

Attach("FrameSymmetry.m");
SetVerbose("FrameSymmetry",0);

function lines(g)
linelist:=[* *];
MPS:=MaxPS(g);
for i:=1 to #MPS do
P,Q,H,N,roots,rho,J,Xi,G:=paraData(g,MPS[i]);
a:=&meet[Eigenspace(gen,1):gen in P].1;
T:=RightTransversal(G,P);
v:=[a^g:g in T];
V:=Parent(v[1]); subs:={@ @}; vecs:=[ ];
for j:=1 to #v do
S:=sub<V|v[j]>;
if S notin subs then
Include(~subs,S);Append(~vecs,v[j]);
end if;
end for;
Append(~linelist,subs);
end for;
return subs;
end function;

function linesi(g,i)
MPS:=MaxPS(g);
P,Q,H,N,roots,rho,J,Xi,G:=paraData(g,MPS[i]);
a:=&meet[Eigenspace(gen,1):gen in P].1;
T:=RightTransversal(G,P);
v:=[a^g:g in T];
V:=Parent(v[1]); subs:={@ @}; vecs:=[ ];
for j:=1 to #v do
S:=sub<V|v[j]>;
if S notin subs then
Include(~subs,S);Append(~vecs,v[j]);
end if;
end for;
return vecs;
end function;

function tLinesi(g,i)
L:=linesi(g,i);
V:=Transpose(Matrix(L));
A:=CanonicalGramian(V);
d:=Nrows(V);

dif:=0; t:=0;
while dif eq 0 do
t:=t+1;
if BaseRing(A) eq Rationals() then
c:=&*[(1+i)/(d+i): i in [0..2*(t-1) by 2]];
else c:=Factorial(d-1)*Factorial(t)/Factorial(d+t-1);
end if;

n:=Ncols(A); sum:=0;
for i:=1 to n do
for j:=1 to n do
sum:=sum + (A[i,j]*ComplexConjugate(A[i,j]))^(t);
end for;
end for;

D:=Diagonal(A); diagsum:=0;
for i:=1 to n do
diagsum:=diagsum + D[i]^t;
end for;
diagsum:=diagsum^2;

dif:=sum-c*diagsum;

end while;
return t-1,A;
end function;

function tLines(g)
tvals:=[];Avals:=[* *];
L:=lines(g);
for i:=1 to #L do
V:=Transpose(Matrix(L[i]));
A:=CanonicalGramian(V);
d:=Nrows(V);

dif:=0; t:=0;
while dif eq 0 do
t:=t+1;
if BaseRing(A) eq Rationals() then
c:=&*[(1+i)/(d+i): i in [0..2*(t-1) by 2]];
else
c:=Factorial(d-1)*Factorial(t)/Factorial(d+t-1);
end if;

n:=Ncols(A); sum:=0;
for i:=1 to n do
for j:=1 to n do
sum:=sum + (A[i,j]*ComplexConjugate(A[i,j]))^(t);
end for;
end for;

D:=Diagonal(A); diagsum:=0;
for i:=1 to n do
diagsum:=diagsum + D[i]^t;
end for;
diagsum:=diagsum^2;

dif:=sum-c*diagsum;

end while;
Append(~tvals,t-1);Append(~Avals,A);
end for;
return tvals,Avals;
end function;

function anglesseq(A)
n:=Ncols(A); ang:=[];
for i in [1..n] do
for j in [1..n] do
Append(~ang,A[i,j]*ComplexConjugate(A[i,j]));
end for;
end for;
return ang;
end function;

function angles(A)
ang:=anglesseq(A);
angs:=Set(ang);
return angs;
end function;

function AngCount(A)
n:=Ncols(A); angseq:=anglesseq(A);angs:=IndexedSet(angseq);count:=ZeroSequence(Rationals(),#angs);
for i in [1..#angs] do
for j in [1..#angseq] do
if angs[i] eq angseq[j] then
count[i]:=count[i]+1;
end if;
end for;
end for;
return angs,count;
end function;


function Angles(A)
ang:=angles(A);
return #ang-1;
end function;

function Anglesi(g,i)
t,A:=tLinesi(g,i);
num:=Angles(A);
return num;
end function;

load "cliff.m";