# # copyright 2003 Baldoni-Silva Vergne De Loera. # ######################################################################### # # WELCOME: # The following software counts the number of integral flows on a network with given # excess values. This is the same as the number of lattice points of a flow polytope: # We call this the KOSTANT PARTITION FUNCTION of a flow polytope. # The input data for general networks are # a) the (1,-1,0)-incidence matrix A of the network. matrix A of dimension n+1xN and rank n,with # columns of the form e_i-e_j. We let m_ij denote the multiplicity of e_i-e_j in the matrix A. # b) an element a=(a_1,...,a_n) with integer coordinates that specifies # the excesses at each node # When the polytope is a trasportation polytope or the complete graph, # this data is not necessary, it is enough to have the parameters of the dimensions. # The formula for the Kostant partitition function that we will compute # is the one given in Theorem 18 of our paper (Counting integer flows on Networks), where # t_i=m_i,i+1 +...+m_i,r+1-1. There are two ingredients in the formula: # (1) The special permutations that contain our vector a and # (2) the computation of certain iterated residues. # HOW TO USE THIS SOFTWARE? Read the examples at the end of the file, they are # very easy to follow and almost self-explanatory. # TO RUN THE EXAMPLES: remove the comment # and read into maple! # # QUESTIONS: write to totalresidue@math.ucdavis.edu # ############################################################################### with(linalg): interface(quiet=true); #_____________________________________________________________________________ #___________________CODE STARTS HERE, EXAMPLES ARE IN LINE 371 _______________ #------------------------------------------------------------------------------ # tipodigrafo contains the information for the network, that is the # computation of the multiplicity, the exponents: graph is for ``generic'' # networks. transpoly is the transportation polytope, completegraph is the # complete graph. # Graph has as input the matrix A that MUST BE GIVEN # AS COORDINATES OF e[i]-e[j] (NO e[n+1]=0 ) so that the number of rows # is the rank-1. For transpoly and complete graph the parameters # giving the dimension are enough without having to write the matrix, # so # The output is always the matrix of the multiplicity, the element # t[i], the number of columns of the matrix, the maximum of the # multiplicity and somma, that is the sum of the roots corresponding to # the matrix columns. #----------------------------------------------------------------------------- graph:=proc(p,q,A) local C,N,n,r,s,st,maxi,sst,su,suu,t,i,qq,d,k,pp,m,tt,f,B,somma,so; sst:=[]; su:=[]; t:=[]; suu:=[]; somma:=[]; N:=coldim(A); r:=rowdim(A); n:=r-1; s:=0; st:=0; for i from 1 to N do s:=add( A[j,i]*e[j], j=1..r ); st:=st+m[op(1,op(1,s)),op(1,op(2,op(2,s)))]; qq:=[op(st)];od; for k from 1 to r do so:=add(A[k,j],j=1..N); somma:=[op(somma),so]; od; for d from 1 to nops(qq) do if op(0,op(d,qq))=m then sst:=[op(sst),[1,[op(1,op(d,qq)),op(2,op(d,qq))]]]; else sst:=[op(sst),[op(1,op(d,qq)),[op(1,op(2,op(d,qq))),op(2,op(2,op(d,qq)))]]] ; fi; od; su:=[seq(op(2,op(i,sst)),i=1..nops(sst))]; suu:=[seq(op(1,op(i,sst)),i=1..nops(sst))];maxi:=max(op(suu)); #print("st:=",st); print("q:=",q); print("su:=",su); print("suu:=",suu); print("sst:=",sst); print(maxi); for k from 1 to n do for pp from 1 to n+1 do if member([k,pp],su,'b') then m[k,pp]:=op(1,op(b,sst)); else m[k,pp]:=0; fi; od; od; for i from 1 to n do tt[i]:=sum(m[i,j],j=i+1..n+1)-1; t:=[op(t),tt[i]];od; f:=(i,j)->m[i,j]; B:=[matrix(n,n+1,f),t,N,maxi,somma]; RETURN(B); end: #----------------------------------------------------------------------- completegraph:=proc(p,q,A) local somma,t,n,m,tt,i,j,f,B; n:=p+q-1;t:=[];tt[0]:=0; for i from 1 to n do for j from i+1 to n+1 do m[i,j]:=1;tt[i]:=(n+1)-i-1; od; for j from 1 to i do m[i,j]:=0; od; t:=[op(t),tt[i]]; od; somma:=[seq(n-2*i+2,i=1..n+1)]; f:=(i,j)->m[i,j];B:=[matrix(n,n+1,f),t,n*(n+1)/2,1,somma]; RETURN(B); end: #---------------------------------------------------------------------------- transpoly:=proc(p,q,A) local somma,n,m,k,t,tt,i,j,f,B;n:=p+q-1; for i from 1 to n do for j from 1 to n+1 do if (member(i,{seq(k,k=1..p)})) and member(j,{seq(k,k=p+1..n+1)}) then m[i,j]:=1; else m[i,j]:=0; fi; od; od; t:=[seq(q-1,d=1..p),seq(-1,d=p+1..n)]; somma:=[seq(1,i=1..p),seq(-1,i=1..q)]; f:=(i,j)->m[i,j]; B:=[matrix(n,n+1,f),t,(p)*(q),1,somma]; RETURN(B); end: #------------------------------------------------------------------------------ # The following procedure checks that the vector a (given by the only first # n coordinates) is indeed inside the cone C(\Delta+), it returns true or false # check_vector:=proc(vector) local c,n,i,s; s:=0; n:=nops(vector); c:=true; i:=1; while i<=n and c=true do s:=s+vector[i]; if s<0 then c:=false; else i:=i+1; fi;od;c; end: #----------------------------------------------------------------------------- # The procedure defvector makes the element a regular if it is not # regular(=inside of a chamber). It is needed to find a chamber # that contain a in its closure. # defvector is by definition a+ epsilon sum(\alpha, \alpha in \Phi) + epsilon^2 (sum e[i]); # the number of coordinates of the vector is p+q-1. defvector:=proc(p,q,A,vector,tipodigrafo) local def,n,ma,m,l ; n:=p+q-1;ma:=tipodigrafo(p,q,A); m:=op(4,ma);l:=op(5,ma); def:=[seq(op(i,vector)+1/(2*(n^2)*m)*l[i]+(1/(2*(n^2)*m))^2,i=1..n)];def; end: #------------------------------------------------------------------------ #computes all the permutations that appear in the iterated residue formula special_permutations:=proc(l,s,vector,tipodigrafo) local perm2,p,q,i,t,set,k,v,w,num,a,b;t:=l+s-1; perm2:=[]; w:=[seq(vector[i],i=1..t)]; #p is the permutation parametrized by [p[1],p[2],p[3],p[4]...] p:=[seq(1,i=1..t)]; # q is the permutation p parametrized the usual way by [q[1],q[2],....] q:=[seq(i,i=1..t)]; set:=[$1..t]; v:=0; i:=1; if check_vector(w)=true then while p[1]<=t do if i=t-1 then a:=set[1]; b:=set[2]; if (tipodigrafo=transpoly) then if ((b<=l) and ((q[i-1]-a)*v<0) and ((a-b)*(v+w[a])<0)) then q[t-1]:=a; q[t]:=b; perm2:=[op(perm2),q]; else fi; if (a<=l )and ((q[i-1]-b)*v<0) and ((b-a)*(v+w[b])<0) then q[t-1]:=b; q[t]:=a; perm2:=[op(perm2),q]; else fi; else if ((q[i-1]-a)*v<0) and ((a-b)*(v+w[a])<0) then q[t-1]:=a; q[t]:=b; perm2:=[op(perm2),q]; else fi; if ((q[i-1]-b)*v<0) and ((b-a)*(v+w[b])<0) then q[t-1]:=b; q[t]:=a; perm2:=[op(perm2),q]; else fi; fi; p[i]:=3; fi; if p[i]>t-i+1 then if i>1 then i:=i-1; k:=q[i]; v:=v-w[k]; set:=[op(1..p[i]-1,set),k,op(p[i]..t-i,set)]; p[i]:=p[i]+1; fi; else k:=op(p[i],set); if (i=1) or ((q[i-1]-k)*v<0) then q[i]:=k; v:=v+w[k]; set:=[op(1..p[i]-1,set), op(p[i]+1..t-i+1,set)]; i:=i+1; p[i]:=1; else p[i]:=p[i]+1; fi; fi; od; fi; perm2; end: #--------------------------------------------------------------------------------------------------- #coeexp computes the coefficients of x^i in the taylor expansion of 1/(1-x)^b.The input are b and i. coeexp:=proc(b,i) local v,j ; v:=1; for j from 1 to i do v:=v*(b+j-1);od;v:=v/i!;end: #-------------------------------------------------------------------------------------------------- #invi computes the taylor expansion up to order u-1 of 1/(1-x)^b, it uses coeexp. invi:=proc(x,b,u) local ss,i,j,c;option remember;for j from 0 to u-1 do c[j]:=coeexp(b,j);od; ss:=(add(c[i]*x^i,i=0..(u-1)));ss;end: #----------------------------------------------------------------------------------------------- # the following procedure computes the Taylor expansion up to order u # of a function constructed by the function g using L. The input are # two integers p,q with p+q-1=rank A ,a the vector defining the # polytope, a permutation L and tipodigrafo that specifies the graph being # considered.If tipodigraph is not graph, then the matrix can be anything. trunc_next_functions:=proc(p,q,A,g,L,u,tipodigrafo,a) local r,C,B,k,n,ff,d,ma,aa; ma:=tipodigrafo(p,q,A);n:=p+q-1; d:=nops(L);k:=L[d]; C:={seq(i,i=1..k-1)} intersect {seq(L[i],i=1..d)}; B:={seq(i,i=(k+1)..n)} intersect {seq(L[i],i=1..d)}; aa[k]:=a[k]+ma[2][k]; ff:=g*invi(-x[k],-aa[k],u)*1/x[k]^(ma[1][k,n+1])*(1/mul(x[j]^ma[1][j,k],j=C))*(1/mul((-x[j])^ma[1][k,j],j=B))*mul(invi(x[k]/x[j],ma[1][j,k],u),j=C)*mul(invi(x[k]/x[j],ma[1][k,j],u),j=B); ff; end: #---------------------------------------------------------------------------------------------- # the procedure E computes an estimate for the order of pole appearing in the residue formula, # independently from the permutation. #$\frac{m(n-k+1)(n-k+2)}{2} -(n-k)$ The imput A can be anything if tipodigrafo is not grafo E:=proc(p,q,A,tipodigrafo) local D,m,s,i,n,ts; s:=[];n:=p+q-1; if tipografo=grafo then D:=grafo(p,q,A); m:=op(4,D); else m:=1; fi; for i from 1 to n do ts:=(m*(n-i+1)*(n-i+2)/2)-(n-i); s:=[op(s),ts]; od; s;end: #----------------------------------------------------------------------- # Eold computes the poles using an old formula. Eold:=proc(p,q) local s,i,n,ts;s:=[];n:=p+q-1; for i from 1 to n do ts:=(n-i+1)*(n-i+2)/2; s:=[op(s),ts]; od; end: # The following procedure computes the same thing as # trunc_next_functions but specify as u the order of the poles # estimates given by procedure E and doing the iterated residue for # the all permutation L. Again if tipodigrafo is not graph, then the # input A can be anything, otherwise it has to be a matrix. RRK:=proc(p,q,A,L,tipodigrafo,a) local EE,tt,n,LL,B,gg,i,K,ma; n:=p+q-1;ma:=tipodigrafo(p,q,A); EE:=E(p,q,A,tipodigrafo);gg:=1;LL:=[op(1..n,L)];i:=1;while i<=n do tt:=trunc_next_functions(p,q,A,gg,LL,EE[n-i+1],ma,a); gg:=sort(expand(coeff(tt,x[LL[n-i+1]],-1))); i:=(i+1);LL:=[op(1..(n-i+1),LL)]; od; gg; end: segnop:=proc(p) local se,i; option remember;se:=1; for i from 1 to nops(p)-1 do if (p[i]>p[i+1]) then se:=se*(-1); fi; od; se; end: # The procedure number_kostant computes the iterated residues for all # the permutations that are in the formula summing it up with the signs # defined by the procedure segno number_kostant:=proc(p,q,A,tipodigrafo,a) local lista,t,n,ma,num,kos,i,g;option remember; n:=p+q-1; lista:=special_permutations(p,q,defvector(p,q,A,a,tipodigrafo),tipodigrafo); ma:=tipodigrafo(p,q,A); kos:=0: num:=nops(lista); i:=1; while i<=num do g:=RRK(p,q,A,lista[i],ma,a);kos:=kos+segnop(lista[i])*g;i:=i+1:od;kos;end: #----------------------------------------------------------- polynomial_kostant:= proc(p,q,A,tipodigrafo,a,v) local lista,t,n,ma,num,kos,i,g;option remember; n:=p+q-1; lista:=special_permutations(p,q,defvector(p,q,A,a,tipodigrafo),tipodigrafo); ma:=tipodigrafo(p,q,A);kos:=0:num:=nops(lista); i:=1; while i<=num do g:=RRK(p,q,A,lista[i],ma,v); kos:=kos+segnop(lista[i])*g; i:=i+1: od; kos; end: #-------------------------EXAMPLES----------------------------------------------------------- # We begin with examples in our paper (see last section in the paper). # EXAMPLE 1 Transportation polytopes (in vector last entry is deleted). #V45:=[3046,5173,6116,10928,-182,-778,-3635,-9558]; # number_kostant(4,5,A,transpoly,V45); # expected answer 23196436596128897574829611531938753 #========================================== # EXAMPLE 2 Another transportation polytope, a big one! V55:=[30201,59791,70017,41731,58270,-81016,-68993,-47000,-43001]; # number_kostant(5,5,A,transpoly,V55); # The answer is 24640538268151981086397018033422264050757251133401758112509495633028, takes time! #special_permutations(5,5,V55, transpoly): nops(S55) #=========================================================== # EXAMPLE 3 Computing the Ehrhart polynomial. # m:=polynomial_kostant(4,0,A,completegraph,[1,2,3],[v1,v2,v3]): factor(m); #================================================================ # EXAMPLE 4 Complete graph K7 #number_kostant(7,0,A,completegraph,[21128,45716,79394,-76028,-31176,66462,-105496]); #polynomial_kostant(7,0,A,completegraph,[21128,45716,79394,-76028,-31176,66462,-105496],[v1,v2,v3,v4,v5,v6,v7]); factor(%); #========================================================== # EXAMPLE 5 A general network #A2:=matrix([[1,0,1,1,0,0],[-1,1,0,0,1,0],[0,-1,-1,0,0,1],[0,0,0,-1,-1,-1]]); #number_kostant1(4,0,A2,graph,[1,2,3]); #10, as in the case complete graph! #m1:=polynomial_kostant(4,0,A2,graph,[1,2,3],[v[1],v[2],v[3]]):factor(m1); #============================================================================== # EXAMPLE 6 Another general network # A5:=matrix([[0,1,1,0,0,0,1,0,1,0],[1,0,0,1,0,1,-1,0,0,0],[-1,-1,0,0,0,0,0,1,0,1],[0,0,-1,-1,1,0,0,0,0,-1],[0,0,0,0,-1,-1,0,-1,-1,0]]); # graph(p,q,A5); # This creates the data structure dealing with that graph. next we print the number of flows: # number_kostant(5,0,A5,graph,[1,2,3,4]); # This happens to be the same as for the complete graph # number_kostant(5,0,A5,completegraph,[1,2,3,4]); # Now we compare the corresponding Ehrhart polynomials. # m4:=polynomial_kostant(4,0,A5,graph,[1,2,3,4],[v[1],v[2],v[3],v[4]]):m5:=factor(m4); # m6:=polynomial_kostant(4,0,A5,completegraph,[1,2,3,4],[v[1],v[2],v[3],v[4]]):m7:=factor(m6);