newPackage( "EDPolytopeCD2", Version => "2.5", Date => "November 7, 2017", Authors => {{Name => "Martin Helmer", Email => "m.helmer@math.ku.dk", HomePage => ""}}, Headline => "Computes ED degrees of toric varieites X_A from the polytope conv(A), B matrix support in codim 2...", DebuggingMode => false, Reload => true ); needsPackage "OldPolyhedra" export{"PolarDegrees", "EDdeg", "CMVolumes", "CMClass", "normalizedVolume", "hyperSimplexVertices", "Output", "TextOutput", "BMatCD2", "faceLatticeAB", "muA", "muB", "degB", "allInrelLine", "indZ", "volB", "PolarDegreesCD2", "EDdegCD2", "CMVolumesCD2", "CMClassCD2" } volB=method(TypicalValue=>ZZ) volB (MutableHashTable):=(alpha)->( Balpha:=alpha#"Balpha"; vec:=first Balpha; primVec:=0; vol:=0; tempDet:=0; if allInrelLine(Balpha) then( primVec=(1/gcd(vec))*vec; ind:=position(vec,i->i!=0); for i from 0 to #Balpha-1 do( tempDet=det(matrix{primVec,Balpha_i}); <<"primVec= "<Boolean) allInrelLine (List) :=(Balpha)->( if #Balpha==1 then return true; --W:=for b in Balpha list 1/(matrix{b}*(transpose matrix{b}))_(0,0)*b; W:=for b in Balpha list (1/gcd(b))*b; v:=#W; --<<"scaled vector= "<ZZ) muA (Matrix,MutableHashTable,MutableHashTable):=(A,alpha,beta)->( return mu(A,alpha,beta); ); muB=method(TypicalValue=>ZZ) muB (List,MutableHashTable,MutableHashTable):=(B,BalphaHash,BbetaHash)->( Balpha:=BalphaHash#"Balpha"; Bbeta:=BbetaHash#"Balpha"; --print "Find mu with"; --<<"Bbeta= "<ZZ) degB (List):=(Balpha)->( B:=matrix Balpha; beta1:=sum select(flatten entries (B_{0}),i->i>0); beta2:=sum select(flatten entries (B_{1}),i->i>0); --<<"beta1= "<ll==wl); Ac=A_L; Bc=B^L; curFace#"Balpha"=entries Bc; Acl=entries transpose Ac; ) else( L2:=for wl in curFace#"Aalpha" list position(Acl,ll->ll==wl); Bcalpha:=Bc^(indsSet-set(L2)); curFace#"Balpha"=entries Bcalpha; ); ifaces=append(ifaces,curFace); ); FacePoset=append(FacePoset,ifaces); ); return FacePoset; ); hyperSimplexVertices=method(TypicalValue=>Matrix) hyperSimplexVertices (ZZ,ZZ):=(d,k)->( hsl:=join(toList(k:1),toList((d-k):0)); return transpose matrix unique permutations(hsl); ); CMClassCD2=method(TypicalValue=>RingElement,Options => {TextOutput=>"Quiet",Output=>RingElement,BMatCD2=>true}) CMClassCD2 (Matrix,QuotientRing):=opts->(A,Chring)->( if opts.Output===HashTable then( return PolarDegrees(A,Chring,Output=>HashTable,TextOutput=>opts.TextOutput,BMatCD2=>opts.BMatCD2); ); return (PolarDegrees(A,Chring,Output=>HashTable,TextOutput=>opts.TextOutput,BMatCD2=>opts.BMatCD2))#"CM class"; ); CMClassCD2 (Matrix):=opts->(A)->( if opts.Output===HashTable then( return PolarDegrees(A,Output=>HashTable,TextOutput=>opts.TextOutput,BMatCD2=>opts.BMatCD2); ); return (PolarDegrees(A,Output=>HashTable,TextOutput=>opts.TextOutput,BMatCD2=>opts.BMatCD2))#"CM class" ); CMClass=method(TypicalValue=>RingElement,Options => {TextOutput=>"Quiet",Output=>RingElement,BMatCD2=>false}) CMClass (Matrix,QuotientRing):=opts->(A,Chring)->( if opts.Output===HashTable then( return PolarDegrees(A,Chring,Output=>HashTable,TextOutput=>opts.TextOutput,BMatCD2=>opts.BMatCD2); ); return (PolarDegrees(A,Chring,Output=>HashTable,TextOutput=>opts.TextOutput,BMatCD2=>opts.BMatCD2))#"CM class"; ); CMClass (Matrix):=opts->(A)->( if opts.Output===HashTable then( return PolarDegrees(A,Output=>HashTable,TextOutput=>opts.TextOutput,BMatCD2=>opts.BMatCD2); ); return (PolarDegrees(A,Output=>HashTable,TextOutput=>opts.TextOutput,BMatCD2=>opts.BMatCD2))#"CM class" ); PolarDegreesCD2=method(TypicalValue=>List,Options => {TextOutput=>"Quiet",Output=>List,BMatCD2=>true}) PolarDegreesCD2 (Matrix):=opts->(A)->( n:=numColumns(A); h:=symbol h; Chring:=ZZ[h]/(h^n); return PolarDegrees(A,Chring,Output=>opts.Output,TextOutput=>opts.TextOutput,BMatCD2=>opts.BMatCD2); ); PolarDegrees=method(TypicalValue=>List,Options => {TextOutput=>"Quiet",Output=>List,BMatCD2=>false}) PolarDegrees (Matrix):=opts->(A)->( n:=numColumns(A); h:=symbol h; Chring:=ZZ[h]/(h^n); return PolarDegrees(A,Chring,Output=>opts.Output,TextOutput=>opts.TextOutput,BMatCD2=>opts.BMatCD2); ); PolarDegrees (Matrix,QuotientRing):=opts->(A,Chring)->( EV:=0; TempEV:=0; FacePoset:={}; if opts.Output===HashTable then( TempEV= CMVolumes(A,TextOutput=>opts.TextOutput,Output=>opts.Output,BMatCD2=>opts.BMatCD2); EV=first TempEV; FacePoset=last TempEV; ) else EV=CMVolumes(A,Output=>opts.Output,TextOutput=>opts.TextOutput,BMatCD2=>opts.BMatCD2); m:=#EV-1; delta:={}; temp:=0; for i from 0 to m do( temp=(-1)^(m-i)*sum(i..m,j->(-1)^(j-i)*binomial(j+1,i+1)*(EV_(m-j))); delta=append(delta, temp); ); EVList:=reverse EV; A=transpose matrix unique entries transpose A; n:=numColumns(A); h:=Chring_0; CMClass:=sum(#(EVList),i->EVList_i*h^(n-1-i)); FirstNonZero:=position (delta,i->(not i==0) ); if opts.Output===List then( <<"The toric variety has degree = "<last(EVList),"dual degree"=>first(delta),"hyperplane class"=>Chring_0,"Chow ring"=>Chring,"ED"=>sum(delta),"polar degrees"=>delta,"CM class"=>CMClass,"CM volumes"=>EVList,"FacePoset"=>FacePoset}; return delHash; ); return delta; ); EDdeg=method(TypicalValue=>ZZ,Options => {TextOutput=>"Quiet",Output=>ZZ,BMatCD2=>false}) EDdeg (Matrix):=opts->(A)->( if opts.Output===HashTable then( return PolarDegrees(A,Output=>opts.Output,TextOutput=>opts.TextOutput,BMatCD2=>opts.BMatCD2); ); return sum PolarDegrees(A,TextOutput=>opts.TextOutput,BMatCD2=>opts.BMatCD2); ); EDdegCD2=method(TypicalValue=>ZZ,Options => {TextOutput=>"Quiet",Output=>ZZ,BMatCD2=>true}) EDdegCD2 (Matrix):=opts->(A)->( if opts.Output===HashTable then( return PolarDegrees(A,Output=>opts.Output,TextOutput=>opts.TextOutput,BMatCD2=>opts.BMatCD2); ); return sum PolarDegrees(A,TextOutput=>opts.TextOutput,BMatCD2=>opts.BMatCD2); ); CMVolumesCD2=method(TypicalValue=>List,Options => {TextOutput=>"Quiet",Output=>List,BMatCD2=>true}) CMVolumesCD2 (Matrix):=opts->(A)->( return CMVolumes(A,Output=>opts.Output, TextOutput=>opts.TextOutput,BMatCD2=>opts.BMatCD2); ); CMVolumes=method(TypicalValue=>List,Options => {TextOutput=>"Quiet",Output=>List,BMatCD2=>false}) CMVolumes (Matrix):=opts->(A)->( AinNum:=numColumns(A); A=transpose matrix unique entries transpose A; if opts.TextOutput!="Quiet" and AinNum!=numColumns(A) then ( print "Warning: Input matrix contains duplicate columns, duplicate columns have been removed"; ); if minors(numRows(A),A)!=ideal(1) then ( if opts.TextOutput!="Quiet" then print "Input matrix does not generate the integer lattice, attempting to build a new matrix which generates the lattice and defines the same toric ideal"; betterA:= transpose gens kernel transpose gens kernel A; if minors(numRows(betterA),betterA)==ideal(1) then( if opts.TextOutput!="Quiet" then print "New matrix generated"; A=transpose hermite transpose betterA; ) else( print "Input matrix does not generate the integer lattice, attempting to build a new matrix which generates the lattice and defines the same toric ideal"; error"Matrix generation failed, please enter a matrix whose columns span the full integer lattice"; return 0; ); ); if rank(A)!=min(numRows(A),numColumns(A)) then ( error "Input matrix expected to have maximal rank"; return 0; ); if rank(A) e == 0) and (all(flatten entries HS2, e -> e <= 0)) then( interiorCols=append(interiorCols,c); ); ); ); --print "done int check"; ); timList=append(timList,first tim2); curFace#"interiorCol"=interiorCols; curFace#"Aalpha"=join(interiorCols,curFace#"verticesList"); if i==0 then curFace#"Aalpha"=Acols; if opts.BMatCD2==true then curFace#"index"=indZ(curFace#"Aalpha"); if opts.BMatCD2==true then( if i==0 then ( L:=for wl in curFace#"Aalpha" list position(entries transpose A,ll->ll==wl); Ac=A_L; Bc=B^L; if opts.TextOutput!="Quiet" then <<"Canonical Forms: Ac= "<ll==wl); Bcalpha:=Bc^(indsSet-set(L2)); curFace#"Balpha"=entries Bcalpha; curFace#"subscripts"=L2; ); ); if i==m then( curVol=1; ) else( if opts.BMatCD2==true then ( if i==0 then( curVol=((m-i)!)*volume(f) ) else( if allInrelLine(curFace#"Balpha") then ( v1:=first curFace#"Balpha"; v1=1/(gcd(v1))*v1; PosSum:=0; for b in entries(B^(curFace#"subscripts")) do( tempdet1:=det(matrix({v1,b})); if tempdet1>0 then PosSum=PosSum+tempdet1; ); curVol=PosSum; ) else curVol=(curFace#"index"); ); ) else( curVol=((m-i)!)*volume(f); ); ); curFace#"volume"=curVol; ifaces=append(ifaces,curFace); ); FacePoset=append(FacePoset,ifaces); ); --<<"time to build face lattice= "<ZZ) normalizedVolume (Polyhedron,ZZ):=(P,n)->( if n>dim(P) then return 0; if n==0 then ( return 1; ) else return ((dim(P))!)*volume(P); ); normalizedVolume (Polyhedron):=(P)->( if dim(P)==0 then return 1 else return ((dim(P))!)*volume(P); ); --------------------------------------------------- --Internal Functions -- --mu=method(TypicalValue=>ZZ) -------------------------------------------------- mu =(A,alpha,beta)->( dBeta:=beta#"dim"; dAlpha:=alpha#"dim"; r:=dAlpha-dBeta; d:=numRows(A); n:=numColumns(A); vbeta:=beta#"verticesMat"; mbeta:=beta#"verticesList"; valpha:=alpha#"verticesMat"; malpha:=alpha#"verticesList"; Atemp:=0; Asort3:=beta#"Aalpha"; Asort2:=toList(set(alpha#"Aalpha")-set(beta#"Aalpha")); Asort1:=toList(set(entries transpose(A))-set(alpha#"Aalpha")); M:=transpose(matrix(join(Asort1,Asort2,Asort3))); if numColumns(A)!=numColumns(M) then ( print "There seems to be an error somewhere"; ); W:=(M); Anew:=transpose hermite(transpose(W)); cs:=n-(#Asort2+#Asort3); ce:=n-#Asort3-1; rowInds:={}; inc:=0; dind:=d-1; for w from 0 to dind do( if flatten(entries((Anew^{dind-w})_(toList((ce+1)..(n-1)))))==toList((#Asort3):0) then( inc=inc+1; rowInds=append(rowInds,dind-w); ); if inc==r then break; ); C:=(Anew_{cs..ce})^rowInds; if rank(C)==0 then return 1; big:=convexHull((transpose(matrix{toList(numRows(C):0)})|C)); C1:=transpose matrix delete(toList(numRows(C):0), entries transpose C); little:=convexHull(C1); vol1:=normalizedVolume(big,r); vol2:=normalizedVolume(little,r); vol:=(vol1-vol2); if numRows(C)!=r then( print "Something may be wrong...we seem to have picked a C matrix with the wrong number of rows"; ); return vol; ); TEST /// {* restart needsPackage "EDPolytope" *} A=transpose(matrix{{0,0,1},{0,7,1},{3,0,1},{5,0,1}}) time pd=EDdeg(A); time pd=PolarDegrees(A,TextOutput=>"All") A=hyperSimplexVertices(4,2); time pd=PolarDegrees(A) assert(pd=={4,12,8,4}); time AHash=PolarDegrees(A,Output=>HashTable); CMVolumes(A); --peek AHash; assert(AHash#"CM volumes"=={12, 12, 8, 4}); assert(AHash#"polar degrees"==pd); assert(AHash#"ED"==28); /// TEST /// {* restart needsPackage "EDPolytope" *} A=transpose matrix{{1,0},{0,2},{3,2},{0,1}}; c1=CMClass(A) assert(c1==CMClass(A,ring c1)); assert(EDdeg(A)==26); Avals=(EDdeg(A,Output=>HashTable)); --peek Avals assert(Avals#"ED"==26); assert(Avals#"dual degree"==7); assert(Avals#"degree"==7); assert(Avals#"polar degrees"=={7,12,7}); assert(Avals#"CM volumes"=={4,9,7}); /// TEST /// {* restart needsPackage "EDPolytope" needsPackage "Polyhedra" *} A=hyperSimplexVertices(4,2) EDdeg A A=matrix{{0,1,1,1,1,0},{0,0,1,1,1,1},{1,1,1,1,0,0}} A=matrix{{0,1,1,1,1,0,0,0,0,0},{0,0,1,1,1,1,0,0,0,0},{1,1,1,1,0,0,1,1,1,1}} A_{0..5} A=matrix{{1,0,1,0,0,0},{0,1,0,1,0,0},{1,1,0,0,0,0},{0,0,1,1,0,0},{0,0,0,0,1,1}} A=matrix{{0,1,1,1,1,0,0,0,0,0},{0,0,1,1,1,1,0,0,0,0},{1,1,1,1,0,0,1,1,1,1}} A=matrix{{0,1,2,1,2,1,2,1,0},{0,0,0,1,1,1,1,1,1},{1,1,1,1,1,1,1,0,0}} A=matrix{{0,1,2,3,4,1,2,3,4,1,2,3,4,1,0},{0,0,0,0,0,1,1,1,1,1,1,1,1,1,1},{1,1,1,1,1,1,1,1,1,1,1,1,1,0,0}} A=matrix{{0,1,2,3,1,2,3,1,2,3,1,0},{0,0,0,0,1,1,1,1,1,1,1,1},{1,1,1,1,1,1,1,1,1,1,0,0}} A=matrix{{0,1,2,3,4,5,1,2,3,4,5,1,2,3,4,5,1,0},{0,0,0,0,0,0,1,1,1,1,1,1,1,1,1,1,1,1},{1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,0,0}} A=transpose matrix unique entries transpose A rank(A) rank(A||matrix{toList(numColumns(A):1)}) A=A||matrix{toList(numColumns(A):1)} A=transpose hermite transpose A F2=faces(1,P) #F2 f2v=for f in F2 list 2*volume(f) sum f2v P=convexHull A B=vertices P faces(1,P) 3!*(volume P) EDdeg A betterA= transpose gens kernel transpose gens kernel A if minors(numRows(betterA),betterA)==ideal(1) then( if opts.TextOutput!="Quiet" then print "New matrix generated"; A=transpose hermite transpose betterA; A=matrix{{0,1,2}} rank A /// TEST /// {* restart needsPackage "EDPolytope" *} A=transpose matrix {{ 2,3,1}, {1 ,2,1}, { 3 ,3,1}, { 4, 4,1},{ 4, 5,1},{1,0,1}} A_{0..4} PolarDegrees A_{1..5} P=convexHull A vertices P contains(P,A_{0}) for f in faces(1,P) list contains(f,A_{2}) F=for f in faces(1,P) list (lift(vertices f,ZZ))|A_{0} C=lift(vertices first faces(1,P),ZZ) PolarDegrees F_2 F_2 I=ToricIdeal A_{1..5} c=codim I kk=ZZ/32749 J=minors(c, jacobian(I)) R=ring I --make conormal var n=numgens(R)-1 degs=join(for i from 0 to n list {1,0},for i from 0 to n list {0,1}) S=kk[gens(ring(I)),v_0..v_n,Degrees=>degs] IS=sub(I,S) K=minors(c+1,matrix{toList(v_0..v_n)}||sub((transpose jacobian I),S)) N=saturate((IS+K),sub(J,S)) --this also is polar degrees multidegree N flatten entries gens sub(ideal(gens(R)),S) N2=ideal(random({0,1},S),random({1,0},S))+N W=eliminate(flatten entries gens sub(ideal(gens()),S),N2) W1=W degree W decompose W /// TEST /// {* restart needsPackage "EDPolytope" needsPackage "ED2Ways" needsPackage "Polyhedra" needsPackage "MLComp" *} A=transpose matrix {{ 2,1,1}, {1 ,2,1}, { 2 ,0,1}, { 1, 1,1},{ 0, 2,1},{1,0,1},{0,1,1}} A_{0..4} minors(3,A) PolarDegrees A P=convexHull A vertices P contains(P,A_{0}) for f in faces(1,P) list contains(f,A_{2}) F=for f in faces(1,P) list PolarDegrees (lift(vertices f,ZZ)) C=lift(vertices first faces(1,P),ZZ) PolarDegrees F_2 F_2 I=ToricIdeal A MLDeg I degree I c=codim I kk=ZZ/32749 J=minors(c, jacobian(I)) codim J numgens ring J R=ring I --make conormal var n=numgens(R)-1 degs=join(for i from 0 to n list {1,0},for i from 0 to n list {0,1}) S=kk[gens(ring(I)),v_0..v_n,Degrees=>degs] IS=sub(I,S) K=minors(c+1,matrix{toList(v_0..v_n)}||sub((transpose jacobian I),S)) N=saturate((IS+K),sub(J,S)) N=saturate((IS+K),sub(ideal product(gens(R)),S)) --this also is polar degrees multidegree N flatten entries gens sub(ideal(gens(R)),S) N2=ideal(random({0,1},S),random({1,0},S))+N W=eliminate(flatten entries gens sub(ideal(gens(R)),S),N) W1=W c=codim W J=minors(c, jacobian(W)) codim J degree W decompose W sub(J,{v_0=>1,v_1=>1,v_2=>1,v_3=>3,v_4=>1,v_5=>1,v_6=>1}) R=QQ[x,y,z] f = x^2*y+x*y^2+x^2+x*y+y^2+x+y codim (ideal(jacobian ideal f)+ideal(f)) factor f MLDeg(I) f = x^2*y+x*y^2+x^2+2*x*y+y^2+x+y factor f degree (ideal(jacobian ideal f)+ideal(f)) I=ToricIdeal(A,{1,1,1,3,1,1,1}) MLDeg1(I) f = x^2*y+x*y^2+x^2+3*x*y+y^2+x+y factor f degree (ideal(jacobian ideal f)+ideal(f)) ///