#-----------------------------------------------------------------
# This maple program computes the BBKLP generating function
# for 4x4 contingency tables of fixed margins.
#
# INPUT: a valid feasible margin vector.
#
# OUTPUT: A rational function.
#-----------------------------------------------------------------

interface(quiet=true):
with(linalg):
with(networks):
with(simplex):

#--------------------------------------------------------
# This procedure Gale dualizes the chambers.
#--------------------------------------------------------

realchambers:=proc(B,cellsofchamber)
local elements,thechamber:

elements:={$1..coldim(B)}:
RETURN(convert(map(x->(elements minus x),cellsofchamber),list)):
end:

#------------------------------------------------------------
# Produces unoriented circuits (just as dependent sets no sign).
#---------------------------------------------------------------

Tounsigned_circuits:=proc(circuits)

RETURN(map(x->op(1,x) union op(2,x),circuits)):

end:

#---------------------------------------------------------------
# Turns the elements of the circuit into two groups of variables
#---------------------------------------------------------------

zcircuits:=proc(circuits)
local nova,c:

nova:=[]:
for c in circuits do
nova:=[op(nova),[map(x->z[x],op(1,c)),map(x->z[x],op(2,c))]]:
od:
RETURN(nova):
end:

#---------------------------------------------------------------------------
# This subroutine finds the coordinate of the vertex of the polytope associated
# with a given base and righthandside of the equation.
#-----------------------------------------------------------------------

coordinates_vertex:=proc(B,righthandside,base)
local const,aux,i,vaux,mataux,y,tempo,tempo2,lista:


aux:=[]: 
for i from 1 to coldim(B) do
  if member(i,base) then
    aux:=[op(aux),y.i]:
  else
    aux:=[op(aux),0]:
  fi:
od:
lista:=convert(base,list):
tempo:=map(x->y.x,lista):
mataux:=submatrix(B,[$1..rowdim(B)],lista):
vaux:=vector(righthandside);
#print(`caca`,linsolve(mataux,vaux),`cama`,lista);
const:=convert(linsolve(mataux,vaux),list);
tempo2:={}:
for i from 1 to nops(const) do
  tempo2:=tempo2 union {op(i,tempo)=op(i,const)}:
od:
RETURN(subs(tempo2,aux));
end:


#-------------------------------------------------------------
# This next subroutine changes a List into a monomial in the
# variables z[i], where the indices run on the size of the list.
#
# EXAMPLE: ListToMonomial(coordinates_vertex(U,[b1,b2,b3],{1,2,6}));
#-------------------------------------------------------------

ListToMonomial:=proc(v)
local tempo,i:
global x:

tempo:=1:
for i from 1 to nops(v) do
tempo:=tempo*(z[i]^op(i,v)):
od:
RETURN(tempo);
end:


#---------------------------------------------------------
# This little subroutine computes the support of a list in fact
#---------------------------------------------------------

support:=proc(l)

local i,neg,pos;
neg:={};
pos:={};
for i from 1 to nops(l) do
  if op(i,l)>0 then
    pos:=pos union {i}:
  elif op(i,l)<0 then
    neg:=neg union {i}:
 fi
od;
RETURN([pos,neg]);
end:

#------------------------------------------------------------
# This procedure goes from the Matrix to the undirected graph
# data structure where we can very fast compute the circuits.
#------------------------------------------------------------

FrommatrixTOgraph:=proc(U)
local edgy,i:
global G:

new(G):
addvertex({$1..rowdim(U)},G):
for i from 1 to coldim(U) do
 edgy:=op(1,support(convert(convert(submatrix(U,[$1..rowdim(U)],[i]),vector),list))):
 addedge(edgy,names=z[i],G):
od:
end:


#----------------------------------------------------------------------
# first an simple auxiliary statement that given a list of indexes that are not
# zero creates a 0,1 vector indicating those indices.

indicator:=proc(c,length) 
local i,v:

v:=[]:
for i from 1 to length do
 if member(i,c) then
  v:=[op(v),1]:
 else
  v:=[op(v),0]:
 fi:
od:
RETURN(v);
end:

indicatordouble:=proc(pos,neg,length) 
local i,v:

v:=[]:
for i from 1 to length do
 if member(i,pos) then
  v:=[op(v),1]:
 elif member(i,neg) then
  v:=[op(v),-1]:
 else
  v:=[op(v),0]:
 fi:
od:
RETURN(v);
end:


#-----------------------------------------------------------------------------
# We compute for a given vertex v of the polytope P_b (with b in the
# chamber) Brions series of its cone K_v (all those polytopes are 
# normally equivalent). This gives one term in Brion's rational function 
# associated to the polytope Pb. The right hand side b is undetermined so far.
# At the same time we are determinining the new adjacent vertices to v!!
#--------------------------------------------------------------------------

BrionsVertexseries:=proc(nodo,circuits,circuitos,circuitsnunu,U)
local d,enotation,auxrhs,numerator, denominator,finalists,e,l,f,adenom,bdenom,bector,vv1,vv2,vvr,entrada,nuevobector:
global G,newneighbors:

base:=op(1,nodo):
bector:=op(2,nodo):
auxrhs:=[]:
for d from 1 to nops(vertices(G))-1 do
auxrhs:=[op(auxrhs),b.d]:
od:

enotation:=map(x->z[x],base):
#print(`esto es enotation`,enotation);
numerator:=ListToMonomial(bector):
denominator:=1:
finalists:={}:
#print(`los circuitos`,circuitos);
for e in (edges(G) minus enotation) do
  # print(`el aristo`,e,`ciclo`, fundcyc(enotation union {e},G));
   if member(fundcyc(enotation union {e},G),circuitos,'l') then
     if member(e,op(1,op(l,circuits))) then
       finalists:=finalists union {op(l,circuitsnunu)}:
     else 
       finalists:=finalists union {[op(2,op(l,circuitsnunu)),op(1,op(l,circuitsnunu))]}:
     fi:
   fi:
od:
for f in finalists do

# This part creates the piece of f for the bb series.


    adenom:=ListToMonomial(indicator(op(1,f),nops(edges(G))));
    bdenom:=ListToMonomial(indicator(op(2,f),nops(edges(G))));
    denominator:=denominator*(1-(adenom/bdenom));

# now we must find the base associated to f

entrada:=min(op(map(x->op(x,bector),op(2,f)))):
#print(`estos son las entradas escojo la mas pequena`,entrada);
vv1:=vector(indicator(op(1,f),nops(edges(G)))):
vv2:=vector(indicator(op(2,f),nops(edges(G)))):
vvr:=matadd(vv1,vv2,entrada,-entrada);
#print(matadd(vvr,vector(bector),1,1),bector,`si si como no`,vvr);
 nuevobector:=convert(matadd(vvr,vector(bector),1,1),list):
 newneighbors:=newneighbors union {[op(1,support(nuevobector)),nuevobector]}:
 od;

 RETURN(numerator/denominator);
end:


FromtopcomTOmaple:=proc(chamber)
local newchamber,c:

newchamber:={}:
for c in chamber do
 newchamber:=newchamber union {map(x->x+1,c)}:
od:
print(nops(newchamber));
RETURN(newchamber):
end:

#----------------------------------------------------------
# This subroutine finds an initial minimal feasible solution which is
# a root for the spanning graph of the enumeration.
#-----------------------------------------------------------

rootenum:=proc(A,b)
local constraints,varvector,value:

varvector:=map(i->x.i,[$1..coldim(A)]):
constraints:=geneqns(A,varvector,vector(b)):
value:=subs(minimize(x1,constraints,NONNEGATIVE),varvector):
print(value);
RETURN([op(1,support(value)),value]):
end:

#------------------------------------------------------------------------------
# This subroutine finds the limit of the univariate BB series when the
# auxiliary parameter t tends to zero. The answer is the number of lattice
# points inside the polytope. This works rightnow only for 4x4 tables as
# the substitution of generic values has to be done for each polytope
# more or less independently.
#------------------------------------------------------------------------------

limiting:=proc(f)
local l,s1,s2,s3,s4,s5,s6,s7,s8,s9,s10,s11,s12,s13,s14,s15,s16,a,cosa:

s1:=subs(z[1]=exp(2*t),f):
s2:=subs(z[2]=exp(23*t),s1):
s3:=subs(z[3]=exp(11*t),s2):
s4:=subs(z[4]=exp(t),s3):
s5:=subs(z[5]=exp(17*t),s4):
s6:=subs(z[6]=exp(3*t),s5):
s7:=subs(z[7]=exp(t),s6):
s8:=subs(z[8]=exp(2*t),s7):
s9:=subs(z[9]=exp(5*t),s8):
s10:=subs(z[10]=exp(7*t),s9):
s11:=subs(z[11]=exp(2*t),s10):
s12:=subs(z[12]=exp(23*t),s11):
s13:=subs(z[13]=exp(11*t),s12):
s14:=subs(z[14]=exp(t),s13):
s15:=subs(z[15]=exp(17*t),s14):
s16:=subs(z[16]=exp(3*t),s15):
print(`I evaluate the rational function now`);
l:=limit(s16,t=0);
print(`the number of contingency tables with given margins`,l);
RETURN(l);
end:

#---------------------------------------------------------------
# Finally we compute the Barvinok-Brion series associated to the
# family of convex polytopes in question as we go along generating 
# the vertices in a DFS fashion.
#---------------------------------------------------------------

PolytopeLatticepointseries:=proc(U,circuits,righthandside)
local i,temp,result,allbases,rhside,V,C,new,allvertices,tobevisited,visited_nodes,nonrepeated:
global G,newneighbors:

# These are just initialization steps.

FrommatrixTOgraph(U):
tobevisited:={rootenum(submatrix(U,[$1..rowdim(U)-1],[$1..coldim(U)]),righthandside)}:
print(tobevisited);
allvertices:={}:
visited_nodes:={}:
V:=submatrix(U,[$1..rowdim(U)-1],[$1..coldim(U)]):
temp:=0:

# Now we run through the general framework of finding series and for each found
# vertex.

while tobevisited <> {} do
 new:=op(1,tobevisited):
 newneighbors:={}:
 temp:=temp+BrionsVertexseries(new,zcircuits(circuits),Tounsigned_circuits(zcircuits(circuits)),circuits,V);
 allvertices:=(allvertices union newneighbors) union {new}:
# print(visited_nodes);
 nonrepeated:=newneighbors minus visited_nodes:
 tobevisited:=(tobevisited minus {new}) union nonrepeated:
 visited_nodes:=visited_nodes union {new}:
 print(`so far I have`,nops(allvertices),`and left`,nops(tobevisited));
od:
   RETURN(limiting(temp));
end:

#------------------SETUP DATA-----------------------

# The matrix U is the incidence matrix of K_33 
U:=matrix([
[1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], 
[0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0], 
[0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0], 
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1], 
[1, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0], 
[0, 1, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0], 
[0,0, 1, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 1, 0], 
[0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 1]]):


circuits:=map(x->{op(1,x),op(3,x)}, 
 [[{2, 11}, {1, 4, 5, 6, 7, 8, 9, 12, 13, 14, 15, 16}, {3, 10}],
[{3, 12}, {1, 2, 5, 6, 7, 8, 9, 10, 13, 14, 15, 16}, {4, 11}], [{2, 12}, {1, 3
, 5, 6, 7, 8, 9, 11, 13, 14, 15, 16}, {4, 10}], [{1, 12}, {2, 3, 5, 6, 7, 8, 
10, 11, 13, 14, 15, 16}, {4, 9}], [{4, 7, 14}, {1, 2, 5, 8, 9, 10, 11, 12, 13,
15}, {3, 6, 16}], [{1, 11}, {2, 4, 5, 6, 7, 8, 10, 12, 13, 14, 15, 16}, {3, 9}
], [{2, 8, 15}, {1, 4, 5, 7, 9, 10, 11, 12, 13, 14}, {3, 6, 16}], [{5, 14}, {1
, 2, 3, 4, 7, 8, 9, 10, 11, 12, 15, 16}, {6, 13}], [{4, 6, 15}, {1, 3, 5, 8, 9
, 10, 11, 12, 13, 14}, {2, 7, 16}], [{1, 6, 11, 16}, {3, 4, 5, 8, 9, 10, 14, 
15}, {2, 7, 12, 13}], [{4, 6, 11, 13}, {2, 3, 5, 7, 9, 12, 14, 16}, {1, 8, 10,
15}], [{3, 9, 14}, {2, 4, 5, 6, 7, 8, 11, 12, 13, 16}, {1, 10, 15}], [{4, 10,
13}, {2, 3, 5, 6, 7, 8, 9, 11, 15, 16}, {1, 12, 14}], [{4, 5, 11, 14}, {1, 3,
6, 7, 10, 12, 13, 16}, {2, 8, 9, 15}], [{3, 5, 10, 16}, {1, 2, 6, 8, 11, 12, 
13, 15}, {4, 7, 9, 14}], [{1, 12, 14}, {3, 4, 5, 6, 7, 8, 10, 11, 13, 15}, {2,
9, 16}], [{3, 8, 9, 14}, {1, 4, 5, 6, 10, 11, 15, 16}, {2, 7, 12, 13}], [{3, 6
, 9, 16}, {1, 2, 5, 8, 11, 12, 14, 15}, {4, 7, 10, 13}], [{3, 6, 12, 13}, {1,
4, 5, 7, 10, 11, 14, 16}, {2, 8, 9, 15}], [{2, 7, 9, 16}, {1, 3, 5, 8, 10, 12,
14, 15}, {4, 6, 11, 13}], [{1, 7, 12, 14}, {3, 4, 6, 8, 9, 10, 13, 15}, {2, 5,
11, 16}], [{1, 6, 12, 15}, {2, 4, 7, 8, 9, 11, 13, 14}, {3, 5, 10, 16}], [{1,
8, 11, 14}, {3, 4, 6, 7, 9, 10, 13, 16}, {2, 5, 12, 15}], [{1, 6, 11, 16}, {2,
3, 7, 8, 9, 12, 13, 14}, {4, 5, 10, 15}], [{1, 12, 15}, {2, 3, 5, 6, 7, 8, 9,
10, 14, 16}, {4, 11, 13}], [{4, 6, 9, 15}, {1, 2, 5, 7, 11, 12, 14, 16}, {3, 8
, 10, 13}], [{1, 7, 10, 16}, {2, 4, 5, 8, 9, 11, 14, 15}, {3, 6, 12, 13}], [{4
, 10, 13}, {1, 3, 5, 6, 7, 8, 11, 12, 14, 15}, {2, 9, 16}], [{2, 7, 9, 16}, {1
, 4, 6, 8, 10, 11, 13, 15}, {3, 5, 12, 14}], [{1, 7, 12, 14}, {2, 4, 5, 6, 9,
11, 15, 16}, {3, 8, 10, 13}], [{4, 11, 14}, {1, 3, 5, 6, 7, 8, 9, 10, 13, 16},
{2, 12, 15}], [{5, 10, 15}, {1, 2, 3, 4, 7, 8, 9, 12, 14, 16}, {6, 11, 13}], [
{1, 8, 10, 15}, {2, 4, 6, 7, 9, 11, 13, 16}, {3, 5, 12, 14}], [{4, 5, 11, 14},
{2, 3, 6, 8, 9, 12, 13, 15}, {1, 7, 10, 16}], [{4, 9, 15}, {1, 2, 5, 6, 7, 8,
10, 11, 14, 16}, {3, 12, 13}], [{1, 11, 16}, {2, 4, 5, 6, 7, 8, 9, 10, 14, 15}
, {3, 12, 13}], [{3, 9, 16}, {1, 2, 5, 6, 7, 8, 10, 12, 14, 15}, {4, 11, 13}],
[{3, 6, 9, 16}, {1, 4, 7, 8, 10, 11, 13, 14}, {2, 5, 12, 15}], [{2, 8, 11, 13}
, {3, 4, 5, 7, 9, 10, 14, 16}, {1, 6, 12, 15}], [{4, 6, 9, 15}, {1, 3, 7, 8, 
10, 12, 13, 14}, {2, 5, 11, 16}], [{1, 10, 16}, {2, 3, 5, 6, 7, 8, 11, 12, 13,
15}, {4, 9, 14}], [{8, 11, 14}, {1, 2, 3, 4, 5, 6, 9, 12, 13, 15}, {7, 10, 16}
], [{4, 6, 9}, {1, 3, 7, 8, 10, 11, 13, 14, 15, 16}, {2, 5, 12}], [{6, 12, 15}
, {1, 2, 3, 4, 5, 8, 9, 11, 13, 14}, {7, 10, 16}], [{11, 16}, {1, 2, 3, 4, 5,
6, 7, 8, 9, 10, 13, 14}, {12, 15}], [{1, 6, 11}, {3, 4, 5, 8, 10, 12, 13, 14,
15, 16}, {2, 7, 9}], [{7, 10, 13}, {1, 2, 3, 4, 5, 8, 11, 12, 14, 16}, {6, 9,
15}], [{4, 10, 15}, {1, 2, 5, 6, 7, 8, 9, 11, 13, 16}, {3, 12, 14}], [{5, 11,
14}, {1, 2, 3, 4, 7, 8, 10, 12, 13, 16}, {6, 9, 15}], [{3, 12, 14}, {1, 4, 5,
6, 7, 8, 9, 10, 13, 15}, {2, 11, 16}], [{5, 10, 15}, {1, 2, 3, 4, 6, 8, 11, 12
, 13, 16}, {7, 9, 14}], [{4, 11, 14}, {1, 2, 5, 6, 7, 8, 9, 12, 13, 15}, {3, 
10, 16}], [{6, 9, 16}, {1, 2, 3, 4, 7, 8, 10, 11, 13, 15}, {5, 12, 14}], [{8,
10, 13}, {1, 2, 3, 4, 6, 7, 9, 11, 15, 16}, {5, 12, 14}], [{2, 12, 15}, {1, 4,
5, 6, 7, 8, 9, 11, 13, 14}, {3, 10, 16}], [{8, 9, 14}, {1, 2, 3, 4, 5, 7, 10,
11, 15, 16}, {6, 12, 13}], [{5, 10, 16}, {1, 2, 3, 4, 7, 8, 9, 11, 14, 15}, {6
, 12, 13}], [{1, 12, 15}, {2, 4, 5, 6, 7, 8, 10, 11, 13, 14}, {3, 9, 16}], [{7
, 10, 13}, {1, 2, 3, 4, 6, 8, 9, 12, 15, 16}, {5, 11, 14}], [{1, 11, 16}, {2,
3, 5, 6, 7, 8, 10, 12, 13, 14}, {4, 9, 15}], [{7, 9, 14}, {1, 2, 3, 4, 5, 8, 
10, 12, 15, 16}, {6, 11, 13}], [{5, 10, 16}, {1, 2, 3, 4, 6, 7, 11, 12, 13, 15
}, {8, 9, 14}], [{3, 5, 10}, {1, 4, 6, 8, 11, 12, 13, 14, 15, 16}, {2, 7, 9}],
[{8, 11, 13}, {1, 2, 3, 4, 6, 7, 9, 10, 14, 16}, {5, 12, 15}], [{8, 9, 15}, {1
, 2, 3, 4, 5, 6, 10, 11, 14, 16}, {7, 12, 13}], [{5, 11, 16}, {1, 2, 3, 4, 6,
7, 10, 12, 13, 14}, {8, 9, 15}], [{8, 11, 14}, {1, 2, 3, 4, 5, 7, 9, 10, 13, 
16}, {6, 12, 15}], [{4, 6, 9}, {2, 3, 5, 7, 11, 12, 13, 14, 15, 16}, {1, 8, 10
}], [{4, 10, 15}, {1, 3, 5, 6, 7, 8, 9, 12, 13, 14}, {2, 11, 16}], [{8, 10, 13
}, {1, 2, 3, 4, 5, 7, 11, 12, 14, 15}, {6, 9, 16}], [{3, 6, 9}, {2, 4, 5, 8, 
11, 12, 13, 14, 15, 16}, {1, 7, 10}], [{1, 7, 10}, {3, 4, 6, 8, 9, 12, 13, 14,
15, 16}, {2, 5, 11}], [{9, 14}, {1, 2, 3, 4, 5, 6, 7, 8, 11, 12, 15, 16}, {10,
13}], [{5, 12, 15}, {1, 2, 3, 4, 6, 8, 10, 11, 13, 14}, {7, 9, 16}], [{1, 6, 
11}, {2, 4, 7, 8, 9, 12, 13, 14, 15, 16}, {3, 5, 10}], [{1, 8, 10}, {3, 4, 6,
7, 9, 11, 13, 14, 15, 16}, {2, 5, 12}], [{8, 10, 15}, {1, 2, 3, 4, 5, 7, 9, 12
, 13, 14}, {6, 11, 16}], [{1, 6, 12}, {2, 3, 7, 8, 9, 11, 13, 14, 15, 16}, {4,
5, 10}], [{5, 11, 16}, {1, 2, 3, 4, 6, 8, 9, 10, 14, 15}, {7, 12, 13}], [{3, 6
, 9}, {1, 4, 7, 8, 10, 12, 13, 14, 15, 16}, {2, 5, 11}], [{8, 11, 13}, {1, 2,
3, 4, 5, 6, 10, 12, 14, 15}, {7, 9, 16}], [{2, 15}, {1, 4, 5, 6, 7, 8, 9, 10,
11, 12, 13, 16}, {3, 14}], [{1, 7, 14}, {3, 4, 6, 8, 9, 10, 11, 12, 13, 16}, {
2, 5, 15}], [{1, 6, 15}, {2, 4, 7, 8, 9, 10, 11, 12, 13, 16}, {3, 5, 14}], [{4
, 6, 13}, {2, 3, 5, 7, 9, 10, 11, 12, 15, 16}, {1, 8, 14}], [{10, 15}, {1, 2,
3, 4, 5, 6, 7, 8, 9, 12, 13, 16}, {11, 14}], [{3, 8, 9}, {2, 4, 5, 6, 10, 11,
13, 14, 15, 16}, {1, 7, 12}], [{4, 7, 9}, {1, 2, 6, 8, 10, 11, 13, 14, 15, 16}
, {3, 5, 12}], [{1, 14}, {3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 15, 16}, {2, 13}],
[{1, 8, 11}, {2, 4, 6, 7, 9, 10, 13, 14, 15, 16}, {3, 5, 12}], [{1, 7, 12}, {2
, 3, 6, 8, 9, 10, 13, 14, 15, 16}, {4, 5, 11}], [{7, 12, 14}, {1, 2, 3, 4, 5,
6, 9, 11, 13, 16}, {8, 10, 15}], [{2, 8, 9}, {1, 3, 6, 7, 11, 12, 13, 14, 15,
16}, {4, 5, 10}], [{9, 16}, {1, 2, 3, 4, 5, 6, 7, 8, 10, 11, 14, 15}, {12, 13}
], [{7, 12, 14}, {1, 2, 3, 4, 5, 8, 9, 10, 13, 15}, {6, 11, 16}], [{11, 13}, {
1, 2, 3, 4, 5, 6, 7, 8, 10, 12, 14, 16}, {9, 15}], [{2, 8, 9}, {3, 4, 5, 7, 10
, 11, 13, 14, 15, 16}, {1, 6, 12}], [{1, 15}, {2, 4, 5, 6, 7, 8, 9, 10, 11, 12
, 14, 16}, {3, 13}], [{2, 8, 11}, {1, 4, 5, 7, 9, 10, 13, 14, 15, 16}, {3, 6,
12}], [{2, 7, 12}, {1, 3, 5, 8, 9, 10, 13, 14, 15, 16}, {4, 6, 11}], [{3, 6, 
13}, {2, 4, 5, 8, 9, 10, 11, 12, 15, 16}, {1, 7, 14}], [{3, 16}, {1, 2, 5, 6,
7, 8, 9, 10, 11, 12, 13, 14}, {4, 15}], [{3, 5, 14}, {1, 4, 6, 8, 9, 10, 11, 
12, 15, 16}, {2, 7, 13}], [{1, 6, 15}, {3, 4, 5, 8, 9, 10, 11, 12, 14, 16}, {2
, 7, 13}], [{2, 8, 11}, {1, 3, 5, 6, 9, 12, 13, 14, 15, 16}, {4, 7, 10}], [{4,
6, 11}, {1, 2, 5, 7, 9, 12, 13, 14, 15, 16}, {3, 8, 10}], [{2, 7, 12}, {1, 4,
5, 6, 9, 11, 13, 14, 15, 16}, {3, 8, 10}], [{1, 16}, {2, 3, 5, 6, 7, 8, 9, 10,
11, 12, 14, 15}, {4, 13}], [{3, 6, 12}, {1, 2, 5, 8, 9, 11, 13, 14, 15, 16}, {
4, 7, 10}], [{10, 16}, {1, 2, 3, 4, 5, 6, 7, 8, 9, 11, 13, 15}, {12, 14}], [{4
, 7, 9}, {2, 3, 5, 6, 10, 12, 13, 14, 15, 16}, {1, 8, 11}], [{3, 8, 9}, {1, 2,
6, 7, 10, 12, 13, 14, 15, 16}, {4, 5, 11}], [{4, 5}, {2, 3, 6, 7, 9, 10, 11, 
12, 13, 14, 15, 16}, {1, 8}], [{4, 6, 13}, {1, 3, 7, 8, 9, 10, 11, 12, 14, 15}
, {2, 5, 16}], [{1, 8, 14}, {3, 4, 6, 7, 9, 10, 11, 12, 13, 15}, {2, 5, 16}],
[{1, 6, 16}, {2, 3, 7, 8, 9, 10, 11, 12, 13, 15}, {4, 5, 14}], [{2, 16}, {1, 3
, 5, 6, 7, 8, 9, 10, 11, 12, 13, 15}, {4, 14}], [{3, 6, 13}, {1, 4, 7, 8, 9, 
10, 11, 12, 14, 16}, {2, 5, 15}], [{1, 7, 16}, {2, 4, 5, 6, 9, 10, 11, 12, 14,
15}, {3, 8, 13}], [{3, 5, 16}, {1, 2, 6, 8, 9, 10, 11, 12, 14, 15}, {4, 7, 13}
], [{2, 8, 13}, {1, 3, 6, 7, 9, 10, 11, 12, 15, 16}, {4, 5, 14}], [{2, 8, 13},
{3, 4, 5, 7, 9, 10, 11, 12, 14, 15}, {1, 6, 16}], [{3, 5}, {2, 4, 6, 8, 9, 10,
11, 12, 13, 14, 15, 16}, {1, 7}], [{8, 15}, {1, 2, 3, 4, 5, 6, 9, 10, 11, 12,
13, 14}, {7, 16}], [{2, 5}, {3, 4, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16}, {1, 6
}], [{8, 13}, {1, 2, 3, 4, 6, 7, 9, 10, 11, 12, 14, 15}, {5, 16}], [{3, 8, 14}
, {1, 4, 5, 6, 9, 10, 11, 12, 13, 15}, {2, 7, 16}], [{5, 15}, {1, 2, 3, 4, 6,
8, 9, 10, 11, 12, 14, 16}, {7, 13}], [{1, 10}, {3, 4, 5, 6, 7, 8, 11, 12, 13,
14, 15, 16}, {2, 9}], [{1, 8, 15}, {2, 4, 6, 7, 9, 10, 11, 12, 13, 14}, {3, 5,
16}], [{1, 7, 16}, {2, 3, 6, 8, 9, 10, 11, 12, 13, 14}, {4, 5, 15}], [{4, 7, 
14}, {1, 3, 5, 6, 9, 10, 11, 12, 13, 16}, {2, 8, 15}], [{4, 6, 15}, {1, 2, 5,
7, 9, 10, 11, 12, 13, 16}, {3, 8, 14}], [{1, 8, 15}, {2, 3, 5, 6, 9, 10, 11, 
12, 14, 16}, {4, 7, 13}], [{4, 5, 15}, {1, 2, 6, 7, 9, 10, 11, 12, 14, 16}, {3
, 8, 13}], [{8, 11}, {1, 2, 3, 4, 5, 6, 9, 10, 13, 14, 15, 16}, {7, 12}], [{6,
15}, {1, 2, 3, 4, 5, 8, 9, 10, 11, 12, 13, 16}, {7, 14}], [{8, 9}, {1, 2, 3, 4
, 6, 7, 10, 11, 13, 14, 15, 16}, {5, 12}], [{5, 11}, {1, 2, 3, 4, 6, 8, 10, 12
, 13, 14, 15, 16}, {7, 9}], [{6, 9}, {1, 2, 3, 4, 7, 8, 11, 12, 13, 14, 15, 16
}, {5, 10}], [{6, 16}, {1, 2, 3, 4, 5, 7, 9, 10, 11, 12, 13, 15}, {8, 14}], [{
4, 7}, {1, 2, 5, 6, 9, 10, 11, 12, 13, 14, 15, 16}, {3, 8}], [{7, 10}, {1, 2,
3, 4, 5, 8, 9, 12, 13, 14, 15, 16}, {6, 11}], [{6, 12}, {1, 2, 3, 4, 5, 7, 9,
11, 13, 14, 15, 16}, {8, 10}], [{4, 6, 9, 15}, {1, 2, 7, 8, 10, 11, 13, 16}, {
3, 5, 12, 14}], [{4, 5, 11, 14}, {1, 3, 6, 8, 9, 10, 15, 16}, {2, 7, 12, 13}],
[{3, 6, 12, 13}, {1, 2, 7, 8, 9, 11, 14, 16}, {4, 5, 10, 15}], [{1, 6, 12, 15}
, {2, 3, 5, 8, 9, 11, 14, 16}, {4, 7, 10, 13}], [{3, 6}, {1, 4, 5, 8, 9, 10, 
11, 12, 13, 14, 15, 16}, {2, 7}], [{4, 6}, {1, 3, 5, 7, 9, 10, 11, 12, 13, 14,
15, 16}, {2, 8}], [{3, 6, 12, 13}, {2, 4, 5, 7, 9, 10, 15, 16}, {1, 8, 11, 14}
], [{4, 5, 11, 14}, {1, 2, 7, 8, 10, 12, 13, 15}, {3, 6, 9, 16}], [{2, 8, 11,
13}, {1, 4, 6, 7, 9, 10, 15, 16}, {3, 5, 12, 14}], [{2, 5, 12, 15}, {1, 3, 6,
8, 10, 11, 13, 16}, {4, 7, 9, 14}], [{4, 6, 11, 13}, {2, 3, 5, 8, 9, 10, 15, 
16}, {1, 7, 12, 14}], [{2, 5, 12, 15}, {1, 4, 6, 7, 9, 11, 14, 16}, {3, 8, 10,
13}], [{2, 5, 11, 16}, {1, 3, 6, 8, 9, 12, 14, 15}, {4, 7, 10, 13}], [{2, 7, 9
, 16}, {3, 4, 5, 6, 10, 12, 13, 15}, {1, 8, 11, 14}], [{1, 6, 11, 16}, {2, 4,
5, 7, 9, 12, 14, 15}, {3, 8, 10, 13}], [{2, 5, 11, 16}, {1, 4, 6, 7, 10, 12, 
13, 15}, {3, 8, 9, 14}], [{4, 6, 11, 13}, {1, 2, 7, 8, 9, 12, 14, 15}, {3, 5,
10, 16}], [{1, 7, 10, 16}, {3, 4, 5, 6, 9, 12, 14, 15}, {2, 8, 11, 13}], [{2,
8, 9, 15}, {1, 4, 6, 7, 11, 12, 13, 14}, {3, 5, 10, 16}], [{1, 8, 10, 15}, {3,
4, 5, 6, 9, 11, 14, 16}, {2, 7, 12, 13}], [{2, 7, 9, 16}, {1, 3, 6, 8, 11, 12,
13, 14}, {4, 5, 10, 15}], [{1, 8, 10, 15}, {2, 3, 5, 6, 11, 12, 13, 16}, {4, 7
, 9, 14}], [{4, 5, 11, 14}, {1, 2, 6, 7, 9, 12, 15, 16}, {3, 8, 10, 13}], [{3,
10, 13}, {2, 4, 5, 6, 7, 8, 9, 12, 15, 16}, {1, 11, 14}], [{2, 8, 9, 15}, {1,
3, 5, 6, 11, 12, 14, 16}, {4, 7, 10, 13}], [{4, 6, 11, 13}, {1, 2, 5, 7, 10, 
12, 15, 16}, {3, 8, 9, 14}], [{3, 6, 9, 16}, {2, 4, 5, 8, 10, 11, 13, 15}, {1,
7, 12, 14}], [{3, 6, 9, 16}, {2, 4, 5, 7, 11, 12, 13, 14}, {1, 8, 10, 15}], [{
2, 8, 9, 15}, {3, 4, 5, 6, 10, 11, 13, 16}, {1, 7, 12, 14}], [{1, 7, 10, 16},
{2, 3, 5, 8, 11, 12, 13, 14}, {4, 6, 9, 15}], [{1, 6, 12, 15}, {2, 4, 5, 7, 10
, 11, 13, 16}, {3, 8, 9, 14}], [{1, 6, 11, 16}, {2, 3, 5, 8, 10, 12, 13, 15},
{4, 7, 9, 14}], [{4, 6, 9, 15}, {2, 3, 5, 7, 10, 12, 13, 16}, {1, 8, 11, 14}],
[{2, 7, 9, 16}, {1, 4, 5, 6, 11, 12, 14, 15}, {3, 8, 10, 13}], [{3, 6, 12, 13}
, {1, 2, 5, 8, 10, 11, 15, 16}, {4, 7, 9, 14}], [{2, 7, 9, 16}, {3, 4, 5, 8, 
10, 11, 13, 14}, {1, 6, 12, 15}], [{2, 8, 11, 13}, {1, 3, 6, 7, 9, 12, 14, 16}
, {4, 5, 10, 15}], [{4, 6, 9, 15}, {1, 3, 5, 8, 10, 11, 14, 16}, {2, 7, 12, 13
}], [{1, 7, 10, 16}, {3, 4, 6, 8, 9, 11, 13, 14}, {2, 5, 12, 15}], [{3, 5, 10,
16}, {1, 4, 6, 8, 9, 11, 14, 15}, {2, 7, 12, 13}], [{1, 7, 10, 16}, {2, 4, 5,
6, 11, 12, 13, 15}, {3, 8, 9, 14}], [{3, 5, 12, 14}, {1, 2, 6, 8, 9, 11, 15, 
16}, {4, 7, 10, 13}], [{1, 10, 16}, {3, 4, 5, 6, 7, 8, 9, 11, 14, 15}, {2, 12,
13}], [{1, 6, 11, 16}, {2, 4, 7, 8, 9, 10, 13, 15}, {3, 5, 12, 14}], [{1, 10,
15}, {3, 4, 5, 6, 7, 8, 9, 12, 14, 16}, {2, 11, 13}], [{2, 5, 11, 16}, {3, 4,
6, 7, 9, 12, 13, 14}, {1, 8, 10, 15}], [{3, 10, 13}, {1, 4, 5, 6, 7, 8, 11, 12
, 14, 16}, {2, 9, 15}], [{4, 5, 11, 14}, {2, 3, 7, 8, 9, 10, 13, 16}, {1, 6, 
12, 15}], [{3, 5, 10, 16}, {2, 4, 6, 7, 9, 12, 13, 15}, {1, 8, 11, 14}], [{1,
7, 12, 14}, {2, 3, 6, 8, 9, 11, 13, 16}, {4, 5, 10, 15}], [{1, 11, 14}, {3, 4,
5, 6, 7, 8, 10, 12, 13, 16}, {2, 9, 15}], [{3, 6, 12, 13}, {1, 4, 7, 8, 9, 10,
14, 15}, {2, 5, 11, 16}], [{3, 9, 14}, {1, 4, 5, 6, 7, 8, 10, 12, 15, 16}, {2,
11, 13}], [{2, 8, 9, 15}, {3, 4, 5, 7, 10, 12, 13, 14}, {1, 6, 11, 16}], [{4,
6, 11, 13}, {1, 3, 7, 8, 9, 10, 14, 16}, {2, 5, 12, 15}], [{3, 6, 9, 16}, {1,
4, 5, 7, 10, 12, 14, 15}, {2, 8, 11, 13}], [{1, 8, 11, 14}, {2, 3, 5, 6, 9, 12
, 15, 16}, {4, 7, 10, 13}], [{4, 5, 10, 15}, {1, 2, 6, 7, 11, 12, 13, 16}, {3,
8, 9, 14}], [{4, 9, 14}, {1, 3, 5, 6, 7, 8, 10, 11, 15, 16}, {2, 12, 13}], [{2
, 8, 11, 13}, {1, 3, 5, 6, 10, 12, 15, 16}, {4, 7, 9, 14}]]):

#-------------EXAMPLE--------------------------------------
# This is a table that has appeared several times in the literature
# For example the survey paper by Diaconis and Gangolli (1995).

a:=[220, 215, 93, 64, 108, 286, 71]:
PolytopeLatticepointseries(U,circuits,a);