// Thomas Booker-Price // This code is licensed under CC BY: // https://creativecommons.org/licenses/by/4.0/ // Some of the code is adapted from code by Jan Grabowski, who has also permitted licensing under CC BY. SetColumns(220); Z:=IntegerRing(); Q:=Rationals(); MRing:=function(n) return MatrixRing(IntegerRing(),n); end function; A3:=function(a,b,c) return MRing(3)![0,-a,-c,a,0,-b,c,b,0]; end function; // ################################################################################################################################ IntBool:=function(tf) if tf ne 0 then return true; else return false; end if; end function; BoolInt:=function(tf) if tf then return 1; else return 0; end if; end function; VarDegree:=function(monomial, varnumber) // find the power of x(varnumber) in given monomial return TotalDegree(monomial)-TotalDegree(Evaluate(monomial, varnumber, 1)); end function; SplitMonomial:=function(monomial, rank) // given x1^a1.x2^a2...xn^an, return [a1, a2, ..., an] return [VarDegree(monomial, i): i in [1..rank] ]; end function; ComputeMonDegree:=function(monomial, grading, rank) // compute the degree of a monomial given a grading V:=VectorSpace(Rationals(),rank); return (V!SplitMonomial(monomial, rank), V!grading); end function; ComputeDegree:=function(clustervar, grading, rank) // compute the degree of a cluster variable given a grading homognum:=Terms(Numerator(clustervar))[1]; denom:=Denominator(clustervar); return ComputeMonDegree(homognum, grading,rank)-ComputeMonDegree(denom, grading, rank); end function; MatrixB:=function(n) // generate Cartan matrix B_n list:=[0:i in [1..n^2]]; list[n^2-1]:=-2; list[2]:=1; if n ge 3 then for i in [2..n-1] do list[n*(i-1)+1+(i-2)]:=-1; list[n*(i-1)+1+i]:=1; end for; end if; return MatrixRing(IntegerRing(),n)!list; end function; MatrixC:=function(n) // generate Cartan matrix C_n B:=MatrixB(n); B[n,n-1]:=-1; B[n-1,n]:=2; return B; end function; BipartiteB:=function(n) // generate bipartite version of Cartan matrix B_n if n lt 3 then return "Bipartite matrix B_n is not defined for this value of n. Use a larged value."; end if; rows:=[[0^^n]^^n]; rows[1][2]:=[1]; for i in [2..n-1] do rows[i][i-1]:=(-1)^(i-1); rows[i][i+1]:=(-1)^(i-1); end for; rows[n][n-1] := (-1)^(n-1)*2; return MatrixRing(IntegerRing(),n)!rows; end function; BipartiteC:=function(n) // generate bipartite version of Cartan matrix C_n if n lt 3 then return "Bipartite matrix C_n is not defined for this value of n. Use a larged value."; end if; rows:=[[0^^n]^^n]; rows[1][2]:=[1]; for i in [2..n-1] do rows[i][i-1]:=(-1)^(i-1); rows[i][i+1]:=(-1)^(i-1); end for; rows[n-1][n]:=rows[n-1][n]*2; rows[n][n-1]:= (-1)^(n-1); return MatrixRing(IntegerRing(),n)!rows; end function; // functions for matrix algebra CA structure ###################################################################################### MatCoordRing:=function(k,l) // generate matrix coordinate ring of rank kxl - adapted from code by Jan Grabowski (adapted from MatrixCA.magma) IndexSet:={ [i,j] : i in [1..k], j in [1..l] }; IndexSetSeq:=Sort(SetToSequence(IndexSet)); IndexSetSeq:=Sort(SetToSequence(IndexSet)); IndexStringSeq:=[]; for S in IndexSetSeq do // get a list of strings the form "11", "12", "13" etc, up to "kl" s:=[]; for i in [1..2] do s[i]:=IntegerToString(S[i]); end for; Append(~IndexStringSeq,&* s); end for; NameSeq:=[]; for S in IndexStringSeq do // list of strings "x_ij" to be used as generator names Append(~NameSeq,"x" cat S); end for; A:= RationalFunctionField(Q,#NameSeq); AssignNames(~A,NameSeq); return A; end function; MACoordMat:=function(k,l) // generate kxl coordinate matrix A:=MatCoordRing(k,l); return Matrix(A,k,l, [ A.i : i in [1..k*l] ] ); end function; VN:=function(i,j, k,l) // returns the number of the MA quiver's vertex in position (i,j), when arranged in a grid - frozen vertices are given the highest numbers if j ne l then return (i-1)*l-(i-1)+j; else return k*l-(i-1); end if; end function; MAExMat:=function(k,l) // generate exchange matrix for matrix coordinate algebra M:=ZeroMatrix(Z,k*l,k*l); for i in [2..(k-1)] do // "interior" vertices for j in [2..(l-1)] do thisV:=VN(i,j,k,l); M[thisV][VN(i-1,j-1,k,l)]:=1; M[thisV][VN(i+1,j,k,l)]:=1; M[thisV][VN(i,j+1,k,l)]:=1; M[thisV][VN(i,j-1,k,l)]:=-1; M[thisV][VN(i-1,j,k,l)]:=-1; M[thisV][VN(i+1,j+1,k,l)]:=-1; end for; end for; for i in [2..(k-1)] do // "left boundary interior" thisV:=VN(i,1,k,l); M[thisV][VN(i+1,1,k,l)]:=1; M[thisV][VN(i,2,k,l)]:=1; M[thisV][VN(i-1,1,k,l)]:=-1; M[thisV][VN(i+1,2,k,l)]:=-1; end for; M[1][VN(1,2,k,l)]:=1; M[1][VN(2,1,k,l)]:=1; M[1][VN(2,2,k,l)]:=-1; // "top left" M[VN(k,1,k,l)][VN(k-1,1,k,l)]:=-1; // "bottom left" for j in [2..(l-1)] do // "lower boundary interior" thisV:=VN(k,j,k,l); M[thisV][VN(k-1,j-1,k,l)]:=1; M[thisV][VN(k-1,j,k,l)]:=-1; end for; M[VN(k,l,k,l)][VN(k-1,l-1,k,l)]:=1; // "bottom right" for j in [2..(l-1)] do // "top interior" thisV:=VN(1,j,k,l); M[thisV][VN(1,j+1,k,l)]:=1; M[thisV][VN(2,j,k,l)]:=1; M[thisV][VN(1,j-1,k,l)]:=-1; M[thisV][VN(2,j+1,k,l)]:=-1; end for; M[VN(1,l,k,l)][VN(1,l-1,k,l)]:=-1; // "top right" for i in [2..(k-1)] do // "right boundary interior" thisV:=VN(i,l,k,l); M[thisV][VN(i-1,l-1,k,l)]:=1; M[thisV][VN(i,l-1,k,l)]:=-1; end for; return M; end function; MARows:=function(r,s,k,l) // row set for minor corresponding to variable at vertex (r,s) of cluster. See Def 5.1 of Jan's Kent notes. return SetToSequence({(k-r+1)..(k-r+s)} meet {1..k}); end function; MACols:=function(r,s,k,l) return SetToSequence({(l-s+1)..(l-s+r)} meet {1..l}); end function; MAInitCl:=function(k,l) // initial cluster for kxl matrix algebra initcluster:=[]; for i in [1..k] do for j in [1..l] do initcluster[VN(i,j,k,l)]:=Determinant(Submatrix( MACoordMat(k,l), MARows(i,j,k,l), MACols(i,j,k,l) )); end for; end for; return initcluster; end function; MAInitGr:=function(k,l) grading:=[]; for i in [1..k] do for j in [1..l] do grading[VN(i,j,k,l)]:=Min(i,j); end for; end for; return grading; end function; // functions for grassmannians CA structure (Convention: new "top left" vertex is given the highest vertex number)############################################### GAExMat:=function(k,l) // Convention: input (k,l) for Gr(k,k+l) M:=ZeroMatrix(Z,k*l+1,k*l+1); MA:=MAExMat(k,l); for i in [1..k*l] do for j in [1..k*l] do M[i][j]:=MA[i][j]; end for; end for; M[k*l+1][1]:=1; M[1][k*l+1]:=-1; return M; end function; GAInitGr:=function(k,l) return [1 : i in [1..k*l+1]]; end function; //################################################################################################################################################ /* Exchange relations (a) Exchange monomials from B [use square B coerced into M! Check with Transpose(B)+B and Transpose(B)*L] */ ExchangeMonomialsSeq:=function(B,r, rank) //returns vectors M1 and M2 i.e. selects B_r^+ and B_r^- and adds the -1's - adapted from code by Jan Grabowski V:=VectorSpace(Rationals(),rank); l:=NumberOfRows(B); M1:=Zero(V); M2:=Zero(V); for i in ({1 .. l} diff {r}) do if (B[i][r] gt 0) then M1[i]:=B[i][r]; M2[i]:=0;// # elif (B[i][r] lt 0) then M1[i]:=0; //# M2[i]:=(-1)*B[i][r]; elif (B[i][r] eq 0) then // # M1[i]:=0;// # M2[i]:=0;// # end if; end for; M1[r]:=-1; M2[r]:=-1; Mons:=[M1,M2]; return Mons; end function; ConvertSeqToMon:=function(cluster,r,seq, rank) // takes a vector (a1,...,an) and a (rational fn) cluster and returns corresponding (rational fn) monomial x^a1...x^an IGNORING THE -1 A:=RationalFunctionField(Rationals(),rank); mon:=Identity(A); for i in [1..(NumberOfColumns(seq))] do if ((Z ! seq[i]) ge 0) then mon:=mon*(A ! cluster[i]^(Z ! seq[i])); end if; end for; return mon; end function; ExchangeMonomials:=function(cluster,B,r, rank) // returns the exchange monomials (rational functions) for mutatation of cluster in direction r given B return [ConvertSeqToMon( cluster,r,ExchangeMonomialsSeq(B,r,rank)[1], rank ),ConvertSeqToMon(cluster,r,ExchangeMonomialsSeq(B,r,rank)[2], rank)]; end function; /* (b) Relation (= coefficients and monomials) */ ExchangeRelation:=function(cluster,B,r) // adapted from code by Jan Grabowski mons:=ExchangeMonomials(cluster,B,r); return [* cluster[r], mons[1]+mons[2] *]; end function; /* Matrix mutations Input: (extended square) exchange matrix B E matrix [only use epsilon=+1, since mutation of B indep. of this] */ EMatrix:=function(B,k, rank) // adapted from code by Jan Grabowski M:=MatrixRing(IntegerRing(),rank); rows:=NumberOfRows(B); E:=[]; for i in [1..rows] do for j in [1..rows] do if not(j eq k) then if (i eq j) then Append(~E,1); else Append(~E,0); end if; elif ((i eq j) and (i eq k)) then Append(~E,-1); else Append(~E,Maximum(0,-1*B[i][k])); end if; end for; end for; return elt; end function; /* F matrix */ FMatrix:=function(B,k, rank) // adapted from code by Jan Grabowski M:=MatrixRing(IntegerRing(),rank); cols:=NumberOfColumns(B); F:=[]; for i in [1..cols] do for j in [1..cols] do if not(i eq k) then if (i eq j) then Append(~F,1); else Append(~F,0); end if; elif ((i eq j) and (i eq k)) then Append(~F,-1); else Append(~F,Maximum(0,B[k][j])); end if; end for; end for; return elt; end function; /* Matrix mutation */ //MatrixMutation:=function(B,k,rank) // return EMatrix(B,k,rank)*B*FMatrix(B,k,rank); //end function; MatrixMutation:=function(B,k,rank:mutable_cols:=rank) // adapted from code by Jan Grabowski mutB:=B; for i in [1..Nrows(B)] do for j in [1..Ncols(B)] do if (i eq k) or (j eq k) then mutB[i][j]:=-B[i][j]; else; if (i le mutable_cols) or (j le mutable_cols) then mutB[i][j]:=B[i][j]+ (Max(B[i][k],0))*B[k][j] + B[i][k]*(Max(-B[k][j],0)); end if; end if; end for; end for; return mutB; end function; MatrixPathMutation:=function(rank,B,path :mutable_cols:=rank) // adapted from code by Jan Grabowski current_matrix:=B; for k in Reverse(path) do current_matrix:=MatrixMutation(current_matrix, k, rank); end for; return current_matrix; end function; /* All cluster variables */ MutatedClusterVariable:=function(cluster,B,k,rank :multi:=false,fast:=false) // cluster is allowed to be a (multi)degree cluster - adapted from code by Jan Grabowski if fast eq false then mons:=ExchangeMonomials(cluster,B,k,rank); m:=(mons[1]+mons[2])/cluster[k]; else mon:=ExchangeMonomialsSeq(B,k,rank)[1]; if multi eq false then m:=&+[cluster[i]*mon[i]: i in [1..rank] ] ; // in the fast version the "clusters" are just degree clusters (remember mon WILL have the -1 entry) else // TODO if needed. Currently this is dealt with in MutatedCluster end if; end if; return m; end function; MutatedCluster:=function(cluster,B,k, rank :multi:=false, fast:=false) // Return mutated cluster in direction k. The variable cluster is allowed to be a multi grading. - adapted from code by Jan Grabowski if multi then grading_vectors:=[]; for i in [1..#cluster] do grading_vectors[i]:=Insert(cluster[i],k,k, [MutatedClusterVariable(cluster[i],B,k,rank: fast:=true)] ); end for; return grading_vectors; else return Insert(cluster,k,k,[MutatedClusterVariable(cluster,B,k, rank :multi:=multi, fast:=fast)]); end if; end function; MutatedClusterAndExchangeMatrix:=function(clusterplusmat,k,rank :multi:=false, fast:=false, mutable_cols:=rank) // adapted from code by Jan Grabowski cluster:=clusterplusmat[1]; B:=clusterplusmat[2]; mutatedclusterplusmat:=[* MutatedCluster(cluster,B,k,rank :multi:=multi, fast:=fast),MatrixMutation(B,k,rank: mutable_cols:=mutable_cols) *]; return mutatedclusterplusmat; end function; IteratedClusterMutation:=function(clusterplusmat,mutationlist, rank :multi:=false, fast:=false, mutable_cols:=rank) //returns [cluster, mat] or [cluster, mat, degrees list] if fast=true - adapted from code by Jan Grabowski l:=#mutationlist; revmutlist:=Reverse(mutationlist); degreeslist:=[]; for i in [1..l] do mutatedclusterplusmat:=MutatedClusterAndExchangeMatrix(clusterplusmat,revmutlist[i], rank :multi:=multi,fast:=fast, mutable_cols:=mutable_cols); clusterplusmat:=mutatedclusterplusmat; if fast eq true and not multi then degreeslist[i]:=mutatedclusterplusmat[1][revmutlist[i]]; end if; end for; if fast eq true then Append(~mutatedclusterplusmat, degreeslist); end if; return mutatedclusterplusmat; end function; NewCVinIteratedClMut:=function(clusterplusmat,mutationlist) // adapted from code by Jan Grabowski return IteratedClusterMutation(clusterplusmat,mutationlist)[1][mutationlist[1]]; end function; Vsum:=function(vector1, vector2) return [vector1[i]+vector2[i]: i in [1..#vector1]]; end function; Vscale:=function(scalar, vector) //scale a "vector" return [scalar*i:i in vector]; end function; ComponentwiseMax:=function(vector1, vector2) return [Max(vector1[i], vector2[i]): i in [1..#vector1]]; end function; DenominatorVectorMutation:=function(rank,dvectors,B,direction) newdvectors:=dvectors; sum1:=[0: i in [1..rank]]; for i in [1..rank] do sum1:= Vsum(sum1, Vscale( Max(B[i,direction],0), dvectors[i] ) ); //[b^t_{ik}]_+ d_{i;t} from (7.7) of CAIV end for; sum2:=[0: i in [1..rank]]; for i in [1..rank] do sum2:= Vsum(sum2, Vscale( Max(-B[i,direction],0), dvectors[i] ) ); //[-b^t_{ik}]_+ d_{i;t} from (7.7) of CAIV end for; newdvectors[direction]:=Vsum( Vscale(-1,dvectors[direction]) , ComponentwiseMax(sum1, sum2)); return newdvectors; end function; e:=function(i,n) vector:=[0:i in [1..n]]; vector[i]:=1; return vector; end function; DvectorPathMutation:=function(rank, B, path) initdvectors:=[Vscale(-1,e(i,rank)) : i in [1..rank]]; currentmatrix:=B; currentdvectors:=initdvectors; for i in [1..#path] do currentdvectors:=DenominatorVectorMutation(rank,currentdvectors,currentmatrix,path[#path-i+1]); currentmatrix:=MatrixMutation(currentmatrix,path[#path-i+1],rank); end for; return currentdvectors; end function; CheckGrading:=function(grading, B: multi:=false, mutable_cols:=Ncols(B)) multigrading:=grading; if multi eq false then multigrading:=[grading]; end if; for j in [1..#multigrading] do checklist:=[ (Vector(multigrading[j]), Transpose(B)[i]): i in [1..Min(#Rows(B), mutable_cols)]]; // UPTO here. Not working "Runtime error in 'Vector': Sequence universe is not a ring" if &or[IntBool(i): i in checklist] then return false; end if; end for; return true; end function; PathMutation:=function(rank,B,path ,grading: multi:=false, fulloutput:=false, compute_dvectors:=true, stringoutput:=true, show_degreecluster:=true, show_matrix:=true, show_heading:=true, show_deglist:=false, custom_initclust:=[], mutable_cols:=rank) // only use slow version if fulloutput is true grading_valid:=CheckGrading(grading,B:multi:=multi,mutable_cols:=mutable_cols); if not grading_valid then return "Error: grading is not homogeneous."; end if; // set up the initial cluster T:=RationalFunctionField(Rationals(),rank); NameSeq:=[]; // TODO change to InitialClusterSetup for i in [1..rank] do Append(~NameSeq,"x" cat IntegerToString(i)); end for; AssignNames(~T,NameSeq); initclust:=[]; if fulloutput eq true then if custom_initclust eq [] then for i in [1..rank] do Append(~initclust,T.i); end for; else initclust:=custom_initclust; end if; else; initclust:=grading; //old: for i in [1..rank] do Append(~initclust,grading[i]); end for; end if; clusterplusmat:=[* initclust, B*]; // carry out the mutations and calculation of degrees degreecl_mat_deglist:=IteratedClusterMutation([* grading, B *], path, rank :multi:=multi,fast:=true, mutable_cols:=mutable_cols); // calculate by mutating degree clusters even if using full variable mutation, to get degrees (fast) if fulloutput eq true then result:=IteratedClusterMutation(clusterplusmat, path, rank :fast:=false, mutable_cols:=mutable_cols); else result:=degreecl_mat_deglist; end if; // only calculate using genuine cluster variables (slow) if fulloutput is true if compute_dvectors eq true then dvectors:=DvectorPathMutation(rank,B,path); else dvectors:=""; end if; // deal with output if fulloutput eq true then print "\n" cat Sprint(result[1]) cat "\n" cat Sprint(result[2]); elif show_matrix eq true then print Sprint(result[2]) cat "\n \n" cat Sprint(dvectors); end if; if show_heading eq true then heading:="Degree \t Mutation Route" cat "\n"; else heading:=""; end if; if fulloutput eq true then output:=Sprint(degreecl_mat_deglist[1][path[1]]) cat "\t" cat Sprint(path); if show_degreecluster eq true then output:= output cat "\n Current degree cluster: \n" cat Sprint(degreecl_mat_deglist[1]); end if; else if not multi then output:=Sprint(result[1][path[1]]) cat "\t" cat Sprint(path); // new degree \t path else new_multidegree:=[result[1][i][path[1]]: i in [1..#result[1]] ]; output:= Sprint(new_multidegree) cat "\t" cat Sprint(path); end if; if show_degreecluster eq true then output:= output cat "\n Current degree cluster: \n" cat Sprint(result[1]); end if; end if; if show_deglist eq true then output:= output cat "\n" cat "Degree list: " cat Sprint(degreecl_mat_deglist[3]); end if; if stringoutput eq true then return heading cat output; else return result; end if; end function; // PathMutation but only need to input k and l, for matrix algebras MAPathMutation:=function(k,l,path: fulloutput:=false, compute_dvectors:=false, stringoutput:=true, show_degreecluster:=true, show_matrix:=true, show_heading:=true, show_deglist:=false) return PathMutation(k*l,MAExMat(k,l), path, MAInitGr(k,l) : fulloutput:=fulloutput, compute_dvectors:=compute_dvectors, stringoutput:=stringoutput, show_degreecluster:=show_degreecluster, show_matrix:=show_matrix, show_heading:=show_heading, show_deglist:=show_deglist, custom_initclust:=MAInitCl(k,l), mutable_cols:=(k*l-l-k+1) ); end function; // PathMutation but only need to input k and l, for grassmannians. Input (k,l) for Gr(k,k+l) GAPathMutation:=function(k,l,path: fulloutput:=false, compute_dvectors:=false, stringoutput:=true, show_degreecluster:=true, show_matrix:=true, show_heading:=true, show_deglist:=false) return PathMutation(k*l+1,GAExMat(k,l), path, GAInitGr(k,l) : fulloutput:=fulloutput, compute_dvectors:=compute_dvectors, stringoutput:=stringoutput, show_degreecluster:=show_degreecluster, show_matrix:=show_matrix, show_heading:=show_heading, show_deglist:=show_deglist, custom_initclust:=MAInitCl(k,l), mutable_cols:=(k*l-l-k+1) ); end function; function RandPath(rank, length: exclude:={}) // generate a random mutation path dirs:={n:n in [1..rank]} diff exclude; path:=[Random(dirs)]; while #path lt length do Append(~path, Random(dirs diff {path[#path]} ) ); end while; return path; end function; function RandPathNoRep(rank, length) // generate a random path in which each mutation direction appears at most once if length gt rank then return "Error: path too long to avoid repitition."; end if; dirs:={n:n in [1..rank]}; path:=[]; while #path lt length do next_dir:=Random(dirs); Append(~path, next_dir ); Exclude(~dirs, next_dir); end while; return path; end function; function RepPath(path, n) return &cat [path : i in [1..n]]; end function; procedure QuickScan(a,b,c) // procedure for getting a very quick idea about a cluster algebra associated to a rank 3 s.s. matrix paths:=[ [1], [2,1], [1,2,1], [2,1,2,1], [1,2,1,2,1], [2,1,2,1,2,1], [1,2,1,2,1,2,1], [2,1,2,1,2,1,2,1], [1,2,1,2,1,2,1,2,1], [2], [1,2], [2,1,2], [1,2,1,2], [2,1,2,1,2], [1,2,1,2,1,2], [2,1,2,1,2,1,2], [1,2,1,2,1,2,1,2], [2,1,2,1,2,1,2,1,2], [3,2], [2,3,2], [3,2,3,2], [2,3,2,3,2], [3,2,3,2,3,2], [2,3,2,3,2,3,2], [3,2,3,2,3,2,3,2], [2,3,2,3,2,3,2,3,2], [2,3], [3,2,3], [2,3,2,3], [3,2,3,2,3], [2,3,2,3,2,3], [3,2,3,2,3,2,3], [2,3,2,3,2,3,2,3], [3,2,3,2,3,2,3,2,3], [3], [1,3], [3,1,3], [1,3,1,3], [3,1,3,1,3], [1,3,1,3,1,3], [3,1,3,1,3,1,3], [1,3,1,3,1,3,1,3], [3,1,3,1,3,1,3,1,3], [3,1], [1,3,1], [3,1,3,1], [1,3,1,3,1], [3,1,3,1,3,1], [1,3,1,3,1,3,1], [3,1,3,1,3,1,3,1], [1,3,1,3,1,3,1,3,1], [1,2,3], [1,2,3,1,2,3], [1,2,3,1,2,3,1,2,3], [1,2,3,1,2,3,1,2,3,1,2,3], [1,2,3,1,2,3,1,2,3,1,2,3,1,2,3], [3,2,1], [3,2,1,3,2,1], [3,2,1,3,2,1,3,2,1], [3,2,1,3,2,1,3,2,1,3,2,1], [3,2,1,3,2,1,3,2,1,3,2,1,3,2,1], [2,1,3], [2,1,3,2,1,3], [2,3,1], [2,3,1,2,3,1], [3,1,2], [3,1,2,3,1,2], [1,3,2], [1,3,2,1,3,2] ]; initexchmat:=A3(a,b,c); grading:=[-b,c,-a]; for i in paths do PathMutation(3,initexchmat, i, grading : show_degreecluster:=false, show_matrix:=false, show_heading:=false); end for; end procedure; /* AllClusterVariables (including initial cluster variables) - adapted from code by Jan Grabowski */ //ExchangeTree:=[* *]; //AllClusterVariablesList:=[* *]; //AllClusterVariablesMutationList:=[* *]; // Example of new usage: AllClusterVariables(3, MatrixB(3), "B": grading:=[2,0,1]); AllClusterVariables:=function(rank,initexchmat,CAtype : grading:=[0..rank]) // adapted from code by Jan Grabowski // things which previously depended on predefined rank ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ A:=RationalFunctionField(Rationals(),rank); NameSeq:=[]; for i in [1..rank] do Append(~NameSeq,"x" cat IntegerToString(i)); // generate a list of names for the indeterminates end for; AssignNames(~A,NameSeq); // assign "x1",...,"xn" to be the names of the indeterminates in A M:=MatrixRing(IntegerRing(),rank); V:=VectorSpace(Rationals(),rank); //x:=function(i) // return A.i; //end function; initclust:=[]; for i in [1..rank] do Append(~initclust,A.i); end for; // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ almostpos:=#(PositiveRoots(IrreducibleRootSystem(CAtype,rank)))+rank; allclvars:={@ A| @}; //use an ordered set to keep things in order vs mutation routes allclmutations:=[* *]; treedepth:=1; tree:=[* *]; t:=Cputime(); for i in [1..rank] do Include(~allclvars,initclust[i]); Append(~allclmutations,[-1*i]); end for; //print [* 0,almostpos,Cputime(t) *]; tree[1]:=[* 1, [* [* initclust,initexchmat,[Z|] *] *] *]; tree[treedepth+1]:=[* 2,[* *] *]; while ((#(allclvars)) lt (almostpos)) do for j in [1..#(tree[treedepth][2])] do //tree[treedepth][2] is list of [cluster, exmat, route] triples at this depth (why not just use j in tree[treedepth][2] ?) for k in [1..rank] do // k will range over directions to mutate in, we have chosen a triple j r:=#(allclvars); ic:=tree[treedepth][2][j][1]; // list of clusters (1st elt of j'th triple in the branch at this depth) ie:=tree[treedepth][2][j][2]; // exchange matrix of (2nd elt of j'th triple in the branch at this depth) im:=tree[treedepth][2][j][3]; // [Z|] (mutation tree route to get here) (3rd elt of j'th triple in the branch at this depth) if ((#im eq 0) or not(k eq im[1])) then // if k wasn't the last direction we mutated in mc:=MutatedCluster(ic,ie,k,rank); me:=MatrixMutation(ie,k,rank); mm:=[k] cat im; for m in [1..rank] do // add to the list of all cluster variables. (should only have to do this once?) Include(~allclvars,mc[m]); end for; Append(~tree[treedepth+1][2],[* mc,me,mm *]); // add this new triple to the next branch if (#(allclvars) eq r+1) then // if there is a new cluster variable Append(~allclmutations,mm); if 2 eq 2 then print [* mc,mm,me*]; end if; end if; if (#allclvars eq almostpos) then // if we have found enough variables, stop break; end if; end if; if (#allclvars eq almostpos) then // break; end if; end for; if (#allclvars eq almostpos) then // break; end if; end for; treedepth:=treedepth+1; tree[treedepth+1]:=[* treedepth+1,[* *] *]; AllClusterVariablesList:=allclvars; AllClusterVariablesMutationList:=allclmutations; //print [* treedepth-1,almostpos-(#AllClusterVariablesList),Cputime(t) *]; if (#allclvars eq almostpos) then // does this do anything? continue; end if; end while; ExchangeTree:=tree; AllClusterVariablesList:=[i: i in allclvars]; AllClusterVariablesMutationList:=allclmutations; if grading eq [0..rank] then AllClusterVarsAndMutations:=[* AllClusterVariablesList,AllClusterVariablesMutationList,Cputime(t) *]; else AllClusterVarsAndMutations:=[* [ : i in AllClusterVariablesList],AllClusterVariablesMutationList,Cputime(t) *]; end if; return AllClusterVarsAndMutations; end function; ACVRadius:=function(rank,initexchmat,radius : grading:=[0..rank]) // adapted from code by Jan Grabowski // things which previously depended on predefined rank ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ A:=RationalFunctionField(Rationals(),rank); NameSeq:=[]; for i in [1..rank] do Append(~NameSeq,"x" cat IntegerToString(i)); // generate a list of names for the indeterminates end for; AssignNames(~A,NameSeq); // assign "x1",...,"xn" to be the names of the indeterminates in A M:=MatrixRing(IntegerRing(),rank); V:=VectorSpace(Rationals(),rank); initclust:=[]; for i in [1..rank] do Append(~initclust,A.i); end for; // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ allclvars:={@ A| @}; //use an ordered set to keep things in order vs mutation routes allclmutations:=[* *]; treedepth:=1; tree:=[* *]; t:=Cputime(); for i in [1..rank] do Include(~allclvars,initclust[i]); Append(~allclmutations,[-1*i]); end for; tree[1]:=[* 1, [* [* initclust,initexchmat,[Z|] *] *] *]; tree[treedepth+1]:=[* 2,[* *] *]; while (treedepth lt radius+1) do for j in [1..#(tree[treedepth][2])] do //tree[treedepth][2] is list of [cluster, exmat, route] triples at this depth for k in [1..rank] do // k will range over directions to mutate in, we have chosen a triple j r:=#(allclvars); ic:=tree[treedepth][2][j][1]; // list of clusters (1st elt of j'th triple in the branch at this depth) ie:=tree[treedepth][2][j][2]; // exchange matrix of (2nd elt of j'th triple in the branch at this depth) im:=tree[treedepth][2][j][3]; // [Z|] (mutation tree route to get here) (3rd elt of j'th triple in the branch at this depth) if ((#im eq 0) or not(k eq im[1])) then // if k wasn't the last direction we mutated in mc:=MutatedCluster(ic,ie,k,rank); me:=MatrixMutation(ie,k,rank); mm:=[k] cat im; for m in [1..rank] do // add to the list of all cluster variables. Include(~allclvars,mc[m]); end for; Append(~tree[treedepth+1][2],[* mc,me,mm *]); // add this new triple to the next branch if (#(allclvars) eq r+1) then // if there is a new cluster variable Append(~allclmutations,mm); //print [* mc,mm,me*]; end if; end if; end for; end for; treedepth:=treedepth+1; tree[treedepth+1]:=[* treedepth+1,[* *] *]; AllClusterVariablesList:=allclvars; AllClusterVariablesMutationList:=allclmutations; end while; AllClusterVariablesList:=[i: i in allclvars]; AllClusterVariablesMutationList:=allclmutations; if grading eq [0..rank] then AllClusterVarsAndMutations:=[* AllClusterVariablesList,AllClusterVariablesMutationList,Cputime(t) *]; else AllClusterVarsAndMutations:=[* [ : i in AllClusterVariablesList],AllClusterVariablesMutationList,Cputime(t) *]; end if; return AllClusterVarsAndMutations; end function; function PrintACV(rank,initexchmat,CAtype : grading:=[0..rank]) // PrintACV(rank,initexchmat,CAtype : grading:=[0..rank]) grad:=grading; if grading ne [0..rank] then checklist:=[ (Vector(grading), Transpose(initexchmat)[i]): i in [1..#Rows(initexchmat)]]; if &or[IntBool(i): i in checklist] then return "Error: grading is not homogeneous."; end if; end if; result:=AllClusterVariables(rank,initexchmat,CAtype : grading:=grad); clvars:=result[1]; clvars:=[i: i in clvars]; mutseqs:=result[2]; if grading eq [0..rank] then heading:="Mutation Route \t Cluster Variable"; output:= &cat [Sprint(mutseqs[i]) cat (&cat [" ":j in [ 1..(5+rank*3-#Sprint(mutseqs[i])) ]]) cat Sprint(clvars[i]) cat "\n" : i in [1..#clvars]]; else heading:="Degree \t Mutation Route \t Cluster Variable"; output:= &cat [Sprint(clvars[i][2]) cat "\t" cat Sprint(mutseqs[i]) cat (&cat [" ":j in [ 1..(5+rank*3-#Sprint(mutseqs[i])) ]]) cat Sprint(clvars[i][1]) cat "\n" : i in [1..#clvars]]; end if; print checklist; return heading cat "\n" cat output; end function; function PrintACVR(rank,initexchmat,radius : grading:=[0..rank]) grad:=grading; if grading ne [0..rank] then checklist:=[ (Vector(grading), Transpose(initexchmat)[i]): i in [1..#Rows(initexchmat)]]; if &or[IntBool(i): i in checklist] then return "Error: grading is not homogeneous."; end if; end if; result:=ACVRadius(rank,initexchmat,radius : grading:=grad); clvars:=result[1]; clvars:=[i: i in clvars]; mutseqs:=result[2]; if grading eq [0..rank] then heading:="Mutation Route \t Cluster Variable"; output:= &cat [Sprint(mutseqs[i]) cat (&cat [" ":j in [ 1..(6+radius*3-#Sprint(mutseqs[i])) ]]) cat Sprint(clvars[i]) cat "\n" : i in [1..#clvars]]; else heading:="Degree \t Mutation Route \t Cluster Variable"; output:= &cat [Sprint(clvars[i][2]) cat "\t" cat Sprint(mutseqs[i]) cat (&cat [" ":j in [ 1..(6+radius*3-#Sprint(mutseqs[i])) ]]) cat Sprint(clvars[i][1]) cat "\n" : i in [1..#clvars]]; end if; return heading cat "\n" cat output; end function; function BriefACVR(rank,initexchmat,radius : grading:=[0..rank]) grad:=grading; if grading ne [0..rank] then checklist:=[ (Vector(grading), Transpose(initexchmat)[i]): i in [1..#Rows(initexchmat)]]; if &or[IntBool(i): i in checklist] then return "Error: grading is not homogeneous."; end if; end if; result:=ACVRadius(rank,initexchmat,radius : grading:=grad); clvars:=result[1]; clvars:=[i: i in clvars]; mutseqs:=result[2]; if grading eq [0..rank] then heading:="Mutation Route"; output:= &cat [Sprint(mutseqs[i]) cat (&cat [" ":j in [ 1..(5+rank*3-#Sprint(mutseqs[i])) ]]) cat "\n" : i in [1..#clvars]]; else heading:="Degree \t Mutation Route"; output:= &cat [Sprint(clvars[i][2]) cat "\t" cat Sprint(mutseqs[i]) cat (&cat [" ":j in [ 1..(6+radius*3-#Sprint(mutseqs[i])) ]]) cat "\n" : i in [1..#clvars]]; end if; return heading cat "\n" cat output; end function; CoxeterNumberOfRootSystem:=function(rootsystem) // adapted from code by Jan Grabowski rank:=Rank(rootsystem); h:=(2*(#PositiveRoots(rootsystem)))/rank; return h; end function; CoxeterNumber:=function(type,rank) R:=IrreducibleRootSystem(type,rank); rank:=Rank(R); h:=(2*(#PositiveRoots(R)))/rank; return h; end function; ExponentsOfRootSystem:=function(type,rank) // adapted from code by Jan Grabowski if type eq "A" then e:=[1..rank]; elif type eq "B" then e:=[]; for i in [1..rank] do e[i]:=2*i-1; end for; elif type eq "C" then e:=[]; for i in [1..rank] do e[i]:=2*i-1; end for; elif type eq "D" then e:=[]; for i in [1..(rank-1)] do e[i]:=2*i-1; end for; e[rank]:=rank-1; elif type eq "E" then if rank eq 6 then e:=[1,4,5,7,8,11]; elif rank eq 7 then e:=[1,5,7,9,11,13,17]; elif rank eq 8 then e:=[1,7,11,13,17,19,23,29]; else e:=[]; end if; elif type eq "F" then if rank eq 4 then e:=[1,5,7,11]; else e:=[]; end if; elif type eq "G" then if rank eq 2 then e:=[1,5]; else e:=[]; end if; elif type eq "H" then if rank eq 3 then e:=[1,5,9]; elif rank eq 4 then e:=[1,11,19,29]; else e:=[]; end if; else e:=[]; end if; return e; end function; ClusterNumber:=function(type,rank) // adapted from code by Jan Grabowski e:=ExponentsOfRootSystem(type,rank); h:=CoxeterNumber(type,rank); l:=[]; for i in [1..rank] do l[i]:=(e[i]+h+1)/(e[i]+1); end for; c:=&*l; return c; end function; /* Extra functions */ IsSkewSymmetricMat:=function(M) // adapted from code by Jan Grabowski return IsZero(M+Transpose(M)); end function; IsSymmetricMat:=function(M) // adapted from code by Jan Grabowski return IsZero(M-Transpose(M)); end function; ModOfSquareMatrix:=function(M) // adapted from code by Jan Grabowski rows:=NumberOfRows(M); matring:=MatrixRing(Z,rows); modM:=[]; for i in [1..rows] do for j in [1..rows] do Append(~modM,Abs(M[i][j])); end for; end for; return elt; end function; QuasiCartanOfSkewSymMat:=function(M) // adapted from code by Jan Grabowski matring:=MatrixRing(Z,NumberOfRows(M)); return (matring ! ModOfSquareMatrix(M))+2*Identity(matring); end function; // ####################################################################################################################################################### //RowToSet:=function(r:multiset:=false, ordered_set:=false) // convert matrix row r to a set of entries // newSet:={}; // if multiset then colset:={* *}; end if; // if ordered_set then colset:={@ @}; end if; // for i in [1..Ncols(r)] do // Include(~newSet, r[1][i]); // end for; // return newSet; //end function; ColSubMat:=function(A,cols) // Return submatrix of A formed by cols=[c_1,...,c_n] return Transpose(Matrix(Transpose(A)[cols])); end function; RowSubMat:=function(A,rows) // Return submatrix of A formed by cols=[c_1,...,c_n] return Matrix(A[rows]); end function; ColToSet:=function(col : multiset:=false, ordered_set:=false) // colset:={}; if multiset then colset:={* *}; end if; if ordered_set then colset:={@ @}; end if; for i in [1..Nrows(col)] do Include(~colset,col[i][1]); end for; return colset; end function; RowToSet:=function(row : multiset:=false, ordered_set:=false) return ColToSet(Transpose(row): multiset:=multiset, ordered_set:=ordered_set); end function; MatToSet:=function(A) // convert matrix A into set of its entries (may need edititing in future if to be used on non-square matrices) rows:=[RowSubMat(A,[i]) : i in [1..Nrows(A)]]; return &join {RowToSet(i): i in rows}; end function; MatMax:=function(A) // return maximum entry of square integer matrix A return Max(MatToSet(A)); end function; // ####################################################################################################################################################### // functions for detecting shortest paths s.t. >=2 appears in the exchange matrix for matrix algebras (2015/2016) ######## IncrementPathAux:=function(path, dirs) // given a path with dirs permissible directions, return "next", "larger", path if #path gt 1 then increment_value:=0; can_increment:=false; while can_increment eq false do // check if it's possible to increase the value of the rightmost entry increment_value:=increment_value+1; if (path[#path]+increment_value ne path[#path-1]) and (path[#path]+increment_value le dirs) then can_increment:=true; end if; if (path[#path]+increment_value gt dirs) then break; end if; end while; if can_increment eq true then return path[1..(#path-1)] cat [path[#path]+increment_value]; else if (path[#path] eq dirs) then // "carry" return $$(path[1..(#path-1)], dirs) cat [1]; else //return $$(path[1..(#path-1)], dirs) cat [path[#path]]; if path[(#path-1)] eq dirs then // entry to right will carry return $$(path[1..(#path-1)], dirs) cat [2]; else return $$(path[1..(#path-1)], dirs) cat [1]; end if; end if; end if; else if path[1] lt dirs then return [path[1]+1]; else return [-1]; end if; end if; end function; SmallestPath:=function(length) // return "smallest" path of given length, i.e. [1,2,1,...,1/2] return [(2-(i mod 2)): i in [1..length]]; end function; IncrementPath:=function(path,dirs) newpath:=IncrementPathAux(path,dirs); if newpath[1] eq -1 then // "overflow" has occured return SmallestPath(#path+1); else return newpath; end if; end function; IncrementPathAlt:=function(path,dirs) // Same as increment path but with the rightmost direction being most significant, rather than leftmost newpath:=IncrementPathAux(Reverse(path),dirs); if newpath[1] eq -1 then // "overflow" has occured return Reverse(SmallestPath(#path+1)); else return Reverse(newpath); end if; end function; PathsList:=function(dirs, maxlength :minlength:=1) // Return all paths with dirs directions that are no longer than given length list:=[]; nextpath:=SmallestPath(minlength); while #nextpath le maxlength do Append(~list, nextpath); nextpath:=IncrementPath(nextpath, dirs); end while; return list; end function; MAFindBigMats:=procedure(k,l,length:minlength:=1,custom_paths:=[], minentry:=2) // Find which paths (no longer than given length) result in a matrix with entries greater than 1 dirs:=k*l-k-l+1; if custom_paths eq [] then paths:=PathsList(dirs, length:minlength:=minlength); else paths:=custom_paths; end if; print "Checking " cat Sprint(#paths) cat " paths. Successful paths found:"; for i in [1..#paths] do mat:=PathMutation(k*l,MAExMat(k,l),paths[i],MAInitGr(k,l) :custom_initclust:=MAInitCl(k,l), mutable_cols:=dirs, stringoutput:=false, compute_dvectors:=false, show_matrix:=false)[2]; if MatMax(mat) ge minentry then print Sprint(paths[i]); end if; if (i mod 500) eq 0 then percent_complete:=Floor(100*Real(Index(paths,paths[i]))/#paths); output:= "Progress: " cat Sprint(percent_complete) cat "%%" cat " Checked " cat Sprint(Index(paths,paths[i])); printf output cat &cat ["\b" : i in [1..(#output-1)]]; end if; end for; end procedure; MAFindBigMatsRand:=procedure(k,l,length,num_paths:minentry:=2) dirs:=k*l-k-l+1; paths:=[]; for i in [1..num_paths] do Append(~paths,RandPath(dirs,length)); end for; print "Checking " cat Sprint(#paths) cat " paths. Successful paths found:"; for i in [1..#paths] do mat:=PathMutation(k*l,MAExMat(k,l),paths[i],MAInitGr(k,l) :custom_initclust:=MAInitCl(k,l), mutable_cols:=dirs, stringoutput:=false, compute_dvectors:=false, show_matrix:=false)[2]; if MatMax(mat) ge minentry then print Sprint(paths[i]); end if; if (i mod 100) eq 0 then percent_complete:=Floor(100*Real(Index(paths,paths[i]))/#paths); output:= "Progress: " cat Sprint(percent_complete) cat "%%" cat " Checked " cat Sprint(Index(paths,paths[i])); printf output cat &cat ["\b" : i in [1..(#output-1)]]; end if; end for; end procedure; // ################################ functions investigating which degrees occur in matrix algebras/grassmannians ############################################################ MACollectDegrees:=function(k,l,num_paths,max_depth:min_depth:=2) // Generate a set containing all degrees found in num_paths random paths of length between min/max depth degree_set:={}; for i in [1..num_paths] do new_depth:=Random(max_depth-min_depth)+min_depth; // length min_depth>1 to max_depth new_path:=RandPath(k*l-k-l+1, new_depth); new_degree:=MAPathMutation(k,l, new_path:show_matrix:=false,stringoutput:=false)[1][1]; //if new_degree lt 1 then printf Sprint(new_path); end if; Include(~degree_set, new_degree); end for; //for i in [1..k*l-k-l+1] do // new_degree:=MAPathMutation(k,l,[i])[1][i]; // Include(~degree_set, new_degree); //end for; return degree_set; end function; GACollectDegrees:=function(k,l,num_paths,max_depth:min_depth:=2) // Generate a set containing all degrees found in num_paths random paths of length between min/max depth degree_set:={}; for i in [1..num_paths] do new_depth:=Random(max_depth-min_depth)+min_depth; // length min_depth>1 to max_depth new_path:=RandPath(k*l-k-l+1, new_depth); new_degree:=GAPathMutation(k,l, new_path:show_matrix:=false,stringoutput:=false)[1][1]; Include(~degree_set, new_degree); end for; return degree_set; end function; // ############################################################################################################################################################################ // Aug 2016 // Functions for finding quiver and degree quiver mutation classes. StandardInitialClusterSetup:=function(rank : return_field:=true) T:=RationalFunctionField(Rationals(),rank); NameSeq:=[]; for i in [1..rank] do Append(~NameSeq,"x" cat IntegerToString(i)); // generate a list of names for the indeterminates end for; AssignNames(~T,NameSeq); // assign "x1",...,"xn" to be the names of the indeterminates in A if return_field eq true then return [* [T.i : i in [1..rank]], T *]; else return [T.i : i in [1..rank]]; end if; end function; PermRowCol:=function(A,perm) // Return simultaneous permutation of rows and cols of matrix A by permutation, input in the form perm=[a_1,...,a_n] P:=PermutationMatrix(IntegerRing(),perm); return P*A*Transpose(P); // P*A permutes rows, A*P^T permutes cols end function; PermSeq:=function(seq, perm) // Return any sequence with entries permuted by perm=[a_1,...,a_n] (not an actual permutation). We can have a_i=a_j. (Avoids using isomorphisms between standard and nonstandard perm group when trying to permute arbitrary sequences.) return [seq[i] : i in perm]; end function; ApplyPerm:=function(seq, perm) // Apply perm to each entry of seq (c.f. PermSeq) return [perm[i]: i in seq]; end function; ApplyPermRep:=function(seq,perm,n) // Apply perm n times new_seq:=seq; for i in [1..n] do new_seq:=ApplyPerm(new_seq,perm); end for; end function; InverseSeqPerm:=function(perm)// Given a permutation (as a sequence), returns the inverse return [Index(perm, i) : i in [1..#perm] ]; end function; PermsFixing:=function(n,fix_list) // Return list of permutation sequences which fix entries in fix_list. (Empty sequence works.) perms:=SetToSequence(Permutations({1..n})); for sigma in perms do if &or[ sigma[i] ne i : i in fix_list] then Exclude(~perms,sigma); end if; end for; return perms; end function; PermsDisjoint:=function(n, disjoint_sets) // Returns permutation sequences *sigma* such that if *set* is in *disjoint_sets* and *i* is in *set* then *i* is sent to something inside *set* under *sigma* perms:=SetToSequence(Permutations({1..n})); for sigma in perms do for set in disjoint_sets do if &or [not (sigma[i] in set) : i in set] then Exclude(~perms,sigma); end if; end for; end for; return perms; end function; PermsPattern:=function(pattern_perm) // Return permutations whose entries agree with the nonzero entries of the sequence pattern_perm n:=#pattern_perm; perms:=SetToSequence(Permutations({1..n})); for sigma in perms do if &or [(pattern_perm[i] ne 0 and sigma[i] ne pattern_perm[i]) : i in [1..n]] then Exclude(~perms,sigma); end if; end for; return perms; end function; PermsPatternSet:=function(pattern_perm_sets) // Return permutations whose entries agree with the sets in the sequence pattern_sets if pattern_perm_sets eq [] then return []; end if; n:=#pattern_perm_sets; perms:=SetToSequence(Permutations({1..n})); for sigma in perms do if &or [not (sigma[i] in pattern_perm_sets[i]) : i in [1..n]] then Exclude(~perms,sigma); end if; end for; return perms; end function; FindIdenticalCols:=function(A,B) colsA:=[ColSubMat(A,[i]) : i in [1..Ncols(A)]]; colsB:=[ColSubMat(B,[i]) : i in [1..Ncols(B)]]; identical_cols:={}; for i in [1..Ncols(B)] do if colsA[i] eq colsB[i] then Include(~identical_cols,i); end if; end for; return identical_cols; end function; Occurences:=function(sequence, a) n:=0; return &+[BoolInt(a eq b): b in sequence]; end function; ColClasses:=function(A) // return the column classes of the matrix A having the same multisets up to sign classes:=[]; remaining_cols:={1..Ncols(A)}; while remaining_cols ne {} do current_col_id:=Min(remaining_cols); current_col_multiset:=ColToSet(ColSubMat(A,[current_col_id]): multiset:=true); neg_current_col_multiset:= ColToSet(-ColSubMat(A,[current_col_id]): multiset:=true); current_class:={current_col_id}; Exclude(~remaining_cols, current_col_id); for i in remaining_cols do match_condition:= ColToSet(ColSubMat(A,[i]): multiset:=true) eq current_col_multiset or ColToSet(ColSubMat(A,[i]): multiset:=true) eq neg_current_col_multiset; if match_condition then Include(~current_class, i); Exclude(~remaining_cols,i); end if; end for; Append(~classes,current_class); end while; return classes; end function; ColClassSignature:=function(A) // Return an invariant of the column classes of A return {*#class: class in ColClasses(A)*}; end function; GenericPermPattern:=function(A,B: col_classesA:=[], col_classesB:=[]) // Return a sequence of sets. The entries of any permutation which could possibly make matrices A and B essentially equivalent must be included in the correspnding sets. (For use with PermsPatternSet.) singletons_covered:={}; identicals:=FindIdenticalCols(A,B); colsA:=[ColSubMat(A,[i]) : i in [1..Ncols(A)]]; colsetsA:=[ColToSet(i:multiset:=true) : i in colsA]; colsB:=[ColSubMat(B,[i]) : i in [1..Ncols(B)]]; colsetsB:=[ColToSet(i:multiset:=true) : i in colsB]; negcolsetsB:=[ColToSet(-i:multiset:=true) : i in colsB]; pattern_sets:=[]; for i in [1..Ncols(B)] do if i in identicals then // if columns are identical, add the singleton and skip to the next column of B Append(~pattern_sets,{i}); Include(~singletons_covered,{i}); continue; end if; current_set:={}; for j in [1..Ncols(A)] do if colsetsA[j] eq colsetsB[i] or colsetsA[j] eq negcolsetsB[i] then Include(~current_set,j); end if; end for; if #current_set eq 1 and current_set in singletons_covered then return []; end if; // The column of A matches only one column of B, but a previous col of A matches only one column of B, so cols of A and B can't be matched. if current_set eq {} then return []; end if; // No columns match Append(~pattern_sets,current_set); if #current_set eq 1 then Include(~singletons_covered,current_set); end if; end for; //"Debug: " cat Sprint(pattern_sets); // Use the singletons to reduce the pattern sets for i in [1..#pattern_sets] do set:=pattern_sets[i]; if #set ge 2 then for singleton in singletons_covered do if (singleton meet set) ne {} then Exclude(~pattern_sets[i], Rep(singleton)); if #pattern_sets[i] eq 1 then if pattern_sets[i] in singletons_covered then return []; end if; // This kind of check would be taken care of in the last few lines, but probably better to check early and save execution Include(~singletons_covered,pattern_sets[i]); end if; if pattern_sets[i] eq {} then return []; end if; // Not enough cols of A having this multiset to match up with cols of B of the same multiset. end if; end for; end if; end for; // Check pattern is not overconstrained for set in pattern_sets do if Occurences(pattern_sets,set) gt #set then return []; end if; end for; return pattern_sets; end function; EqualColsCheck:=function(A,B) // Given matrices A,B of same dimension, return [a_1,...,a_n], where a_i != 0 means i'th column set of B equals (a_i)'th column set of A colsA:=[ColSubMat(A,[i]) : i in [1..Ncols(A)]]; colsetsA:=[ColToSet(i:multiset:=true) : i in colsA]; colsB:=[ColSubMat(B,[i]) : i in [1..Ncols(B)]]; colsetsB:=[ColToSet(i:multiset:=true) : i in colsB]; negcolsetsB:=[ColToSet(-i:multiset:=true) : i in colsB]; checklist:=[0: i in [1..Ncols(A)]]; directions_covered:={}; found_matching_col:=false; for i in [1..#colsetsB] do skip_multiset_check:=false; // put a negative entry if colums are actually identical, i.e. for parts of the quiver which have not changed at all if colsA[i] eq colsB[i] then // note i=j here if not (i in directions_covered) then checklist[i]:=-i; Include(~directions_covered, i); skip_multiset_check:=true; found_matching_col:=true; end if; end if; if not skip_multiset_check then for j in [1..#colsetsA] do if not (j in directions_covered) then if (colsetsA[j] eq colsetsB[i]) or (colsetsA[j] eq negcolsetsB[i]) then // this avoids repeated directions in checklist checklist[i]:=j; Include(~directions_covered, j); found_matching_col:=true; break; end if; end if; end for; if not found_matching_col then return checklist; end if; end if; end for; return checklist; end function; EsseqMat:=function(A,B :colclass_signatureA:={* *}, colclass_signatureB:={* *}, display_message:=true) // Given matrices A and B with the same number of columns, return whether they are essentially equivalent and gives sigma s.t. B=sigma(A). if colclass_signatureA ne {* *} and colclass_signatureA ne colclass_signatureB then //"Debug: early negative detection"; return [*false, []*]; end if; if Ncols(A) ne Ncols(B) then return "Error: A and B have a different number of columns."; end if; rank:=Ncols(A); pattern_perm_sets:=GenericPermPattern(A,B); //Gives sigma s.t. B=sigma(A). Alternatively (B,A) gives sigma s.t. sigma(B)=A when applied as simultaneous row/col perm - need to change col_classes correspondingly. //"Debug: pattern perm" cat Sprint(pattern_perms_sets); if pattern_perm_sets eq [] then if display_message then print "A and B are not equivalent."; end if; //"Debug: Not equivalent due to different column sets."; return [*false, []*]; end if; perms:=PermsPatternSet(pattern_perm_sets); //"Debug: number of permutations is " cat Sprint(#perms); if #perms eq 0 then Sprint(pattern_perm_sets) cat "\n A \n" cat Sprint(A) cat "\n B \n" cat Sprint(B); end if; // Apply the permutations to test for equivalance for sigma in perms do sigma_A:=PermRowCol(A,sigma); are_esseq:= B eq sigma_A or -B eq sigma_A; if are_esseq then if display_message then "A and B are essentially equivalent by " cat Sprint(sigma) cat "."; end if; return [* true, sigma *]; end if; end for; //"Debug: not equivalent after exhausting all possible permutations."; return [*false, []*]; end function; NextMutPath:=function(rank, path, closed) // Determines the order of vertex enumeration in FindDegreeQuivers. The variable closed indicates whether we have just found an equivalent degree seed or not. Returns [rank+1] if the input path is closed and there are no more paths left. nextpath:=path; if path eq [] then return [1]; end if; if closed then done:=false; while done eq false do if #nextpath eq 1 then nextpath[1]:= nextpath[1]+1; done := true; else lastdir:=nextpath[1]; nextpath:=nextpath[2..#nextpath]; if lastdir lt rank and lastdir+1 ne nextpath[1] then nextpath:=[lastdir+1] cat nextpath; done:=true; elif lastdir+1 eq nextpath[1] and lastdir+2 le rank then nextpath:=[lastdir+2] cat nextpath; done:=true; end if; end if; end while; else if nextpath[1] eq 1 then nextpath:= [2] cat nextpath; else nextpath:= [1] cat nextpath; end if; end if; return nextpath; end function; FindQuiverClass_old:=function(A) t:=Cputime(); found:=1; rank:=Ncols(A); quiver_list:=[A]; colclass_signature_list:=[ColClassSignature(A)]; current_path:=[]; depth:=0; quiver_depth_list:=[A]; // Note the i-th entry of this sequence corresponds to depth i-1 closed:=false; while NextMutPath(rank, current_path, closed) ne [rank+1] do current_path:=NextMutPath(rank, current_path, closed); Sprint(current_path); next_dir:=current_path[1]; depth:=#current_path; B:=MatrixMutation(quiver_depth_list[depth],next_dir,rank); closed:=false; for i in [1..found] do if EsseqMat(quiver_list[i], B : colclass_signatureA:=colclass_signature_list[i], colclass_signatureB:=ColClassSignature(B), display_message:=false)[1] or EsseqMat(quiver_list[i], -B : colclass_signatureA:=colclass_signature_list[i], colclass_signatureB:=ColClassSignature(B), display_message:=false)[1] then closed:=true; break; end if; end for; if closed eq false then // found new Append(~quiver_list,B); Append(~colclass_signature_list,ColClassSignature(B)); quiver_depth_list[depth+1]:=B; found:=found+1; "Found " cat Sprint(found); end if; end while; Cputime(t); return #quiver_list; end function; UnderlyingSimpleMat:=function(Q) // Given (quiver) matrix Q, return matrix corresponding to the underlying simple graph, i.e. with all nonzero entries replaced by 1 for i in [1..Nrows(Q)] do for j in [1..Ncols(Q)] do if Q[i,j] ne 0 then Q[i,j]:=1; end if; end for; end for; return Q; end function; UnderlyingSimpleGraph:=function(Q) // Given (quiver) matrix Q, return underlying simple graph return Graph< Ncols(Q) | UnderlyingSimpleMat(Q) >; end function; Quiver:=function(Q) // Given (quiver) matrix Q, return the quiver (represented as a simple graph with labelled vertices encoding the edge weights and directions - this is the only way we can use IsIsomorphic, and hence the nauty package, properly) rank:=Nrows(Q); num_layers:=Floor(Log(2,MatMax(Q))) + 1; // create the skeleton quiver_graph:=Digraph; // Initially start with the three vertices of Layer 0 at the bottom. AssignLabels(~quiver_graph, [VertexSet(quiver_graph)!v: v in [1..rank]], [("Layer 0")^^rank] ); for i in [1..num_layers-1] do AddVertices(~quiver_graph, rank); AssignLabels(~quiver_graph, [ VertexSet(quiver_graph)!v: v in [( i*rank+1 )..( i*rank+rank )] ] , [("Layer " cat Sprint(i))^^rank] ); end for; for l in [2..num_layers] do AddEdges(~quiver_graph, { [j + rank*(l-1), j + rank*(l-1)-rank] : j in [1..rank] } ); AddEdges(~quiver_graph, { [j + rank*(l-1)-rank, j + rank*(l-1)] : j in [1..rank] } ); end for; // encode each weighted edge i---w-->j by adding edges so that w is represented in binary up when moving up layers between the corresponding vertices of the skeleton. for i in [1..rank] do current_row:=(RowSubMat(Q,[i])); for j in [1..rank] do current_entry:=current_row[1][j]; // [1] required to turn row into tuple //"Debug: "Sprint(current_entry); if current_entry gt 0 then binary_seq:=Intseq(current_entry,2); for d in [1..#binary_seq] do if binary_seq[d] eq 1 then quiver_graph +:= [i+rank*(d-1), j+rank*(d-1)]; end if; end for; end if; end for; end for; return quiver_graph; end function; FindQuiverClass:=function(A) t:=Cputime(); found:=1; rank:=Ncols(A); quiver_list:=[A]; current_path:=[]; depth:=0; quiver_depth_list:=[A]; // Note the i-th entry of this sequence corresponds to depth i-1. Used to store matrices along the tree to allow for one mutation per new vertex (rather than computing the whole mutation path starting from the initial matrix). closed:=false; while NextMutPath(rank, current_path, closed) ne [rank+1] do current_path:=NextMutPath(rank, current_path, closed); Sprint(current_path); next_dir:=current_path[1]; depth:=#current_path; B:=MatrixMutation(quiver_depth_list[depth],next_dir,rank); closed:=false; for i in [1..found] do if IsIsomorphic(Quiver(quiver_list[i]), Quiver(B) ) then closed:=true; break; elif IsIsomorphic(Quiver(quiver_list[i]), Quiver(-B) ) then closed:=true; break; end if; end for; if closed eq false then // found new Append(~quiver_list,B); quiver_depth_list[depth+1]:=B; found:=found+1; "Found " cat Sprint(found); end if; end while; Cputime(t); return #quiver_list; end function; NautyIsomToPerm:=function(isom, rank, quiver) inv_perm:= [ Index( Vertices(quiver),isom(i) ): i in [1..rank] ]; // Note Index is required, rather than just isom(i), as isom(i) is a graph vertex and can't be coerced into an integer return InverseSeqPerm(inv_perm); // Inverse is needed because isom is an isomorphism between two quivers, not a permutation of a single quiver. end function; PermGrading:=function(perm, grading) return [grading[i] : i in perm]; end function; IsEsseqDegreeQuiver:=function(A, gradA, B, gradB: allow_neg_grad:=false, only_allow_neg_grad:=false, custom_isom:=false, given_isom:=0, multi:=false) // A and B are matrices of quivers, gradA/B are multigradings (seqs of basis vector seqs) unless multi is false if multi eq false then gradA:=[gradA]; gradB:=[gradB]; end if; grading_dimension:=#gradA; rank:=Ncols(A); QA:=Quiver(A); QB:=Quiver(B); QB_op:=Quiver(-B); if not custom_isom then AB_isom,isom:=IsIsomorphic(QA,QB); if not AB_isom then AB_isom,isom:=IsIsomorphic(QA,QB_op); if not AB_isom then return false,_; end if; end if; perm:=NautyIsomToPerm(isom, rank, QA); // perm is s.t. perm(A)=B (up to sign) else perm:=NautyIsomToPerm(given_isom, rank, QA); end if; result:=true; if not allow_neg_grad and not only_allow_neg_grad then for i in [1..grading_dimension] do if not (PermGrading(perm, gradA[i]) eq gradB[i]) then result:=false; break; end if; end for; elif allow_neg_grad and not only_allow_neg_grad then for i in [1..grading_dimension] do if not( (PermGrading(perm, gradA[i]) eq gradB[i]) or (PermGrading(perm, [-z: z in gradA[i]]) eq gradB[i])) then result:=false; break; end if; end for; elif only_allow_neg_grad then for i in [1..grading_dimension] do if not(PermGrading(perm, [-z: z in gradA[i]]) eq gradB[i]) then result:=false; break; end if; end for; end if; return result, perm; end function; FindDegreeQuivers:=function(A,grading : multi:=false) // Find all degree seeds up to essential equivalence. Unless grading space is 1-dimnl, the multigrading is to be given in the form of a sequence of basis vector sequences if multi eq false then grading:=[grading]; end if; grading_dimension:=#grading; rank:=Ncols(A); // Check grading is admissible for j in [1..grading_dimension] do checklist:=[ (Vector(grading[j]), Transpose(A)[i]): i in [1..rank]]; if &or[IntBool(i): i in checklist] then return "Error: grading is not homogeneous."; end if; end for; // Main function t:=Cputime(); found:=1; quiver_list:=[A]; gradings_list:=[grading]; degrees_found:={}; for i in [1..rank] do // Add multidegrees of initial grading in degrees_found. this_degree:=[0^^grading_dimension]; for j in [1..grading_dimension] do this_degree[j]:=grading[j][i]; end for; Include(~degrees_found, this_degree); end for; current_path:=[]; depth:=0; quiver_depth_list:=[A]; // Note the i-th entry of this sequence corresponds to depth i-1. Used to store matrices along the tree to allow one mutation per new vertex (rather than computing the whole mutation path starting from the initial matrix) gradings_depth_list:=[grading]; closed:=false; while NextMutPath(rank, current_path, closed) ne [rank+1] do current_path:=NextMutPath(rank, current_path, closed); next_dir:=current_path[1]; depth:=#current_path; output:="Current direction: " cat Sprint(next_dir) cat " Current depth: " cat Sprint(depth); printf output cat &cat ["\b" : i in [1..(#output)]]; B:=MatrixMutation(quiver_depth_list[depth],next_dir,rank); current_grading:=[]; for i in [1..grading_dimension] do current_grading[i]:=MutatedCluster(gradings_depth_list[depth][i], quiver_depth_list[depth], next_dir, rank : fast:=true); end for; Include(~degrees_found, [ current_grading[i][next_dir] : i in [1..grading_dimension] ] ); closed:=false; for i in [1..found] do found_isom, isom := IsIsomorphic(Quiver(quiver_list[i]), Quiver(B) ); if found_isom then if IsEsseqDegreeQuiver(quiver_list[i], gradings_list[i], (B), current_grading : custom_isom:=true, given_isom:=isom , multi:=true ) then closed:=true; break; end if; else found_isom, isom :=IsIsomorphic(Quiver(quiver_list[i]), Quiver(-B) ); if found_isom then if IsEsseqDegreeQuiver( quiver_list[i], gradings_list[i], (-B), current_grading : custom_isom:=true, given_isom:=isom , multi:=true ) then closed:=true; break; end if; end if; end if; end for; if closed eq false then // found new quiver Append(~quiver_list,B); quiver_depth_list[depth+1]:=B; Append(~gradings_list,current_grading); gradings_depth_list[depth+1]:=current_grading; found:=found+1; "\nFound " cat Sprint(found) cat " degree quivers and the following degrees:"; Sprint(degrees_found); end if; end while; "Degree seeds found: " cat Sprint(#quiver_list); "Found the following set of degrees: " cat Sprint(degrees_found); "Time taken: "; Cputime(t); return [*quiver_list,gradings_list*]; end function; // ################################################################################################################ // Some file operations FILE_WriteSeq:=procedure(filename, seq :overwrite:=false) // Create a text file containing seq (one element per line) PrintFile(filename, seq[1] , "Magma": Overwrite:=overwrite); if #seq ge 2 then for i in seq[2..#seq] do PrintFile(filename, i, "Magma"); // This does add a newline character. (Note: this may not be interpreted on every text editor, but it is present.) Will append if the file already exists. end for; end if; end procedure; FILE_WriteElt:=procedure(filename, element :overwrite:=false) // Create a text file containing elt (should have only one line) PrintFile(filename, element, "Magma" : Overwrite:=overwrite); end procedure; FILE_Append:=procedure(filename, element) // Append an element to a file PrintFile(filename, element, "Magma"); end procedure; FILE_LoadSeq:=function(filename) // Load and return a sequence from an existing file (expects one element of the sequence per line) seq:=[]; file:=Open(filename, "r"); current_line:= Gets(file); while not IsEof(current_line) do current_element:=eval current_line; Append(~seq,current_element); current_line:=Gets(file); end while; return seq; end function; FILE_LoadList:=function(filename) // Load and return a list from an existing file (expects one element of the list per line) list:=[* *]; file:=Open(filename, "r"); current_line:= Gets(file); while not IsEof(current_line) do current_element:=eval current_line; Append(~list,current_element); current_line:=Gets(file); end while; return list; end function; FILE_LoadObject:=function(filename) // Load and return an object from a file containing only that object contents:=Read(filename); object:= eval contents; return object; end function; FILE_LoadGradingsSeq:=function(filename, grad_dim) // Load a sequence of multigradings of dimension grad_dim (each written as one line per basis vector, plus a line at the start and end), for the variables gradings_list and gradings_depth_list in FDQ_RestoreSession seq:=[]; file:=Open(filename, "r"); current_line_partial:=[]; while true do for i in [1..grad_dim+2] do current_line_partial[i]:= Gets(file); // Will be "[ PowerSequence(IntegerRing()) |" when i=1, "[1,0,1,...]" when 1