Macaulay 2 Code

-- Macaulay2 function for computing the covering number, that is, the
-- minimum number of monomials of degree d required to cover the monomials of
-- degree d + 1.
--
-- v1.0
--
-- Ported to Macaulay2 by Ben Babcock 
-- Based on original CoCoA algorithm by E. Carlini, H. T. Ha, and A. Van Tuyl.
------------------------------------------------------------------------------

coveringNumber = (n,d) -> (
if n == 1 then return 1;
if d == 1 then return n;

k := binomial(d + n - 1, n - 1);
S := QQ[x_1..x_n];
T := QQ[z_1..z_(k - n)];

Ad := flatten entries gens((monomialIdeal(gens S))^d);
AdPlus1 := flatten entries gens((monomialIdeal(gens S))^(d + 1));

Powers := for i from 1 to n list x_i^d;
Powers2 := {};
for i from 1 to n do (
for j from 0 to (#Powers - 1) do (
Powers2 = append(Powers2, x_i*Powers#j);
);
);

AdPrime := select(Ad, i -> not member(i,Powers));
AdPlus1Prime := select(AdPlus1, i -> not member(i, Powers2));

W := for i from 0 to (#AdPrime - 1) list (Powers | drop(AdPrime,{i,i}));

Igens := {};
for i from 0 to (#AdPlus1Prime - 1) do (
m := AdPlus1Prime#i;
mQuotients := for j from 1 to n list (if gcd(m, x_j) == 1 then continue; m//x_j);
Igens = append(Igens,
product apply(mQuotients, m -> (
product for i from 0 to (#W - 1) list (if member(m, W#i) then continue; z_(i+1))
))
);
);

ISimpComp := monomialIdeal Igens;
return(k - dim(T/ISimpComp));
);

- Macaulay2 function for computing the spreading number
-- v1.0
--
-- Ported to Macaulay2 by Ben Babcock 
-- Based on original CoCoA algorithm by E. Carlini, H. T. Ha, and A. Van Tuyl.
------------------------------------------------------------------------------

spreadingNumber = (n,d) -> (
if n == 1 or d == 1 then return 1;
if d == 2 then return n;

k := binomial(d + n - 1, n - 1);
S := QQ[x_1..x_n];
T := QQ[z_1..z_k];

-- In the ring S, first determine all monomials of degree d. Second,
-- determine which pair of monomials have a LCM of degree d+1
monoS := (monomialIdeal(gens S))^d;
Sd := flatten entries gens monoS;
idealGens := {};
for i from 1 to k do (
for j from (i + 1) to k do (
m := lcm(Sd#(i - 1), Sd#(j - 1));
degm := first degree m;
if degm == d + 1 then idealGens = append(idealGens, z_i*z_j);
);
);

-- Construct the Stanley-Reisner ring and compute its dimension
ISimpComp := monomialIdeal idealGens;
return(dim(T/ISimpComp));
)

-- Macaulay2 algorithm for finding a greedy lower bound for the spreading number
-- v0.5
-- Ben Babcock 
------------------------------------------------------------------------------

needsPackage "EdgeIdeals";

getMaxIndepSet = method();

getMaxIndepSet = (n,d) -> (
if n == 1 or d == 1 then return 1;
if d == 2 then return n;

k := binomial(d + n - 1, n - 1);
x := local x;
z := local z;

S := QQ[x_1..x_n];
T := QQ[z_1..z_k];

-- In the ring S, first determine all monomials of degree D. Second,
-- determine which pair of monomials have a LCM of degree D+1

monoS := (monomialIdeal(gens S))^d;
Sd := flatten entries gens monoS;
Pairs := {};
for i from 1 to k do (
for j from (i + 1) to k do (
m := lcm(Sd#(i - 1), Sd#(j - 1));
degm := first degree m;
if degm == d + 1 then Pairs = append(Pairs,{z_i, z_j});
);
);

G := graph(T,Pairs);
L := vertices G;
W := {};

while(#L > 0) do (
v := first L;
if #W == 0 then v = L#(random(0, #L - 1))
else (
-- Select a vertex of minimum degree
dV := apply(L, i -> degreeVertex(G, i));
v = L#(position(dV, i -> i == min dV));
);

-- Add that vertex to the set
W = append(W,v);

-- Remove v and its neighbours from L
L = toList(set(L) - ({v} | neighbors(G,v)));
);

return #W;
);


----------------------------------------------------------------------------
----------------------------------------------------------------------------

-- Macaulay2 function for computing the spreading number of S_n(d)

-- This algorithm uses a recursive random greedy algorithm to generate
-- maximal independent sets, which serve as lower bounds. It then checks
-- square-free monomials of the next degree to see if they encode a larger
-- independent set, using the symmetry of S_n(d) to reduce the number
-- of monomials that must be checked.

-- v0.7 (now with less MemoryLeak, I promise)
-- Ben Babcock [bababcoc@lakeheadu.ca]
------------------------------------------------------------------------------

needsPackage "EdgeIdeals";

-- Given a list of vertices, build an independent set
-- Used in the greedyMaxIndepSet algorithm
buildIndepSet = method();
buildIndepSet(Graph, List, List) := (G,L,I) -> (
    v := first L;

    if #I == 0 then v = L#(random(0, #L - 1))
    else (
-- From the remaining vertices, select the vertex with the smallest degree.
-- This way, we are choosing the vertex that will have the "lowest impact"
-- on the number of vertices removed from contention.
for i from 1 when i < #L do (
if degreeVertex(G, L#i) < degreeVertex(G, v) then v = L#i;
);
    );
    I = append(I, v);
    
    -- Remove any vertices adjacent to v from L
    L = toList(set(L) - ({v} | neighbors(G,v)));

    -- Are there any vertices left?
    if #L == 0 then return I;

    return buildIndepSet(G,L,I);
);

----------------------------------------------------------------------------
----------------------------------------------------------------------------

-- Iterate over monomials of each degree until we have found the maximum
-- independent set
--
-- Required Parameters
-- S The ring QQ[x_1..x_n]
-- G The graph S_n(d)
-- initialBound An initial lower bound on the spreading number
-- W The initial maximal independent set
-- Legend A hash table encoding vertices z_i => monomials in S
-- inverseLegend A hash table encoding monomials in S => vertices z_i
--
checkMonomials = (S, G, initialBound, W, Legend, inverseLegend) -> (
    T := ring G;
    use T;

    -- We will use the best bound we have, which is either the initial bound
    -- or the size of the initial set, whatever is largest
    activeDegree := if initialBound > #W then initialBound else (#W + 1);

    -- Get the symmetric group of degree n
    sGroup := symmetricGroup S;

    -- Create the initial generators of the ideal we'll use to test a candidate set
    -- We'll start with the edge ideal and any permutations of our initial set
    genList := permuteMonomials(S, G, sGroup, Legend, inverseLegend, W) | (edges G);

    while true do (
<< endl << "W is " << W << endl;
<< "Active degree is " << activeDegree << endl;

-- Now we iterate through the squarefree monomials of degree activeDegree.
combo := toList(1..activeDegree);
time passed := while true do (
-- Check if the last combination was empty (i.e., no more squarefree
-- monomials of this degree).
if #combo == 0 then break {};

-- Create the corresponding set of vertices from the combination
m := apply(combo, i -> z_i);

-- Check if this set is independent (i.e., does it contain any smaller
-- maximal independent set? If no, does it contain any edges?
if not any(genList, V -> isSubset(V, m)) then break m;

combo = nextCombination(numgens T, combo);
);

if #passed == 0 then break
else (
W = greedyMaxIndepSet(G, passed);
genList = permuteMonomials(S, G, sGroup, Legend, inverseLegend, W) | genList;
activeDegree = #W + 1;
);
    );
    << endl << "W is " << W << endl;
    return W;
);


-----------------------------------------------------------------------------
-----------------------------------------------------------------------------

-- A greedy algorithm that takes a graph and an independent set on that graph as input.
-- The independent set can be empty, in which case the graph will choose a random vertex
-- from which to begin building the set. However, passing a non-empty, non-independent
-- set will result in an error.
--
-- The function outputs a maximal independent set
------------------------------------------------------------------------------

greedyMaxIndepSet = method();
greedyMaxIndepSet(Graph, List) := (G,W) -> (
    use ring G;

    if #W == 0 then W = buildIndepSet(G, vertices G, {})
    else if not isSubset(W, vertices G) then error "Expected second argument to be a list of vertices.";
    m := product W;
    
    I := edgeIdeal G;
    M := monomialIdeal m;
    J := I : M;

    -- If J is the unit ideal, the set isn't independent!
    if J == 1 then error "Expected an independent set.";

    -- If J is generated by monomials of degree 1, we may be done.
    possiblyMaximal := false;
    if max(apply(degrees J, i->i#0)) == 1 then possiblyMaximal = true;

    -- First we find all vertices not in the neighbourhood of W. If this list has
    -- one element, we add it to the set before saying we are done.
    newVertices := toList(set(select(vertices G, i -> degreeVertex(G, i) > 0)) - (W | neighbors(G, W)));

    if possiblyMaximal and #newVertices == 0 then return W;
    if possiblyMaximal and #newVertices == 1 then return(W | newVertices);

    -- Otherwise, we induce a subgraph on the new list of vertices.
    H := inducedGraph(G, newVertices, OriginalRing => true);
    
    -- Select a random vertex of H and build a new independent set
    V = buildIndepSet(H, newVertices, {});
    
    return greedyMaxIndepSet(H, W | V);
);


----------------------------------------------------------------------------
----------------------------------------------------------------------------

-- Generate the next r-combination from a list m of r elements (lexicographical order)
-- m is a subset of {1,..n}
nextCombination = (n,m) -> (
    r := #m;
    if r > n then return {};

    -- Check if the subset is equal to {n - r + 1,...,n}
    if m == toList((n - r + 1)..n) then return {};

    m = new MutableList from m;
    i := r;

    while m#(i - 1) == (n - r + i) do i = i - 1;
    m#(i - 1) = m#(i - 1) + 1;

    for j from (i + 1) to r do m#(j - 1) = m#(i - 1) + j - i;

    return (new List from m);
);


-----------------------------------------------------------------------------
-----------------------------------------------------------------------------

-- Given a list of vertices, apply each permutation from the symmetric group to
-- every vertex. This results in a list of permuted sets of vertices, or rotations and
-- reflections of the set. This operation preserves independence of sets.
permuteMonomials = (S, G, sGroup, Legend, inverseLegend, W) -> (
    L = {};
    for i when i < #sGroup do (
newGens := {};
sigma := sGroup#i;
for j when j < #W do (
-- Use the legend to look up the value of this vertex in its
-- original ring. Factor it.
f := factor Legend#(W#j);
m := 1;

for k when k < #f do (
m = m*(x_(sigma#(index(f#k#0)) + 1)^(f#k#1));
);

-- Translate this monomial back to the graph ring
newGens = append(newGens, inverseLegend#m);
);
L = append(L, newGens);
    );
    return(unique L);
);


----------------------------------------------------------------------
----------------------------------------------------------------------

recursiveTransform = (dMinusTwo, nMinusOne) -> (
    -- dMinusTwo has monomials from S_n(d - 2)
    -- We must adjust them by multiplying by x_1^2
    use ring first dMinusTwo;
    dMinusTwo = apply(dMinusTwo, i -> i*x_1^2);

    -- nMinusOne has monomials from S_{n - 1}(d)
    -- We must adjust the indices of their indeterminates by +1
    use ring first nMinusOne;
    nMinusOne = apply(nMinusOne, i -> (
apply(factor i, f -> (
x_(index(f#0) + 2)^(f#1)
))
    ));
    return (dMinusTwo | nMinusOne);
);


----------------------------------------------------------------------------
----------------------------------------------------------------------------

-- This is the main event.
-- Required Parameters
-- n The number of indeterminates, x_1..x_n
-- d The degree of the monomials that interest us
-- Options
-- InitialSet An initial independent set to use as the seed.
-- Passing a non-independent set results in an error.
-- Default empty.
-- OutputSet Boolean. When true, output the maximum independent
-- set. When false, output just the number. Default false.
spreadingNumber = method(Options => {InitialSet => {}, OutputSet => false});

spreadingNumber(ZZ,ZZ) := opts -> (n,d) -> (
    if n == 1 or d == 1 then return 1;
    if d == 2 then return n;

    k := binomial(d + n - 1, n - 1);
    S := QQ[x_1..x_n];
    T := QQ[z_1..z_k];

    -- In the ring S, first determine all monomials of degree d. Second,
    -- determine which pair of monomials have a LCM of degree d + 1
    use S;
    Sd := flatten entries gens((monomialIdeal(gens S))^d);
    Pairs := {};
    legend := {};
    for i from 1 to k do (
legend = append(legend, {z_i, Sd#(i - 1)});

for j from (i + 1) to k do (
m := lcm(Sd#(i - 1), Sd#(j-1));
degm := first degree m;
if degm == d + 1 then Pairs = append(Pairs,{z_i,z_j});
);
    );

    -- This table will match up generators of T with monomials of S_n(d)
    legend = hashTable legend;
    inverseLegend := applyPairs(legend, (p,q) -> (q,p));

    -- Construct the graph S_n(d)
    G := graph(T, Pairs);

    -- This is the best lower bound we have so far
    initialBound := k//n;

    -- Use a greedy algorithm to construct an initial maximal independent set
    initialSet := apply(opts.InitialSet, i -> inverseLegend#(sub(i,S)));
    initialSet = greedyMaxIndepSet(G, initialSet);

    -- call to new function here (checkMonomials()?)
    finalSet := checkMonomials(S, G, initialBound, initialSet, legend, inverseLegend);
    
    if opts.OutputSet then return apply(finalSet, v -> legend#v)
    else return #finalSet;
);


------------------------------------------------------------------------------
------------------------------------------------------------------------------

-- Return a list of elements of the symmetric group of degree n
symmetricGroup = method();
symmetricGroup(Ring) := R -> (
    -- This line from the Monomial Ideals chapter, by Serkan Hosten and Gregory Smith,
    -- of Computations in Algebraic Geometry with Macaulay2
    L := terms det matrix table(numgens R, gens R, (i,r) -> r^(i + 1));
    sGroup := {};

    -- Construct a hashtable to represent each permutation (starting from 0 instead of 1)
    for i when i < #L do (
elem := {};
f := factor L#i;
for j when j < #f do (
if index f#j#0 =!= null then elem = append(elem, {index f#j#0, f#j#1 - 1});
);
sGroup = append(sGroup, hashTable elem);
    );

    return sGroup;
);


------------------------------------------------------------------------------
------------------------------------------------------------------------------

-- Macaulay2 algorithm for finding an upper bound for the spreading number
-- 
-- Ben Babcock [bababcoc@lakeheadu.ca]
------------------------------------------------------------------------------

needsPackage "EdgeIdeals";
spreadingNumberBound = (n,d) -> (
if n == 1 or d == 1 then return 1;
if d == 2 then return n;

k := binomial(d + n - 1, n - 1);
x := local x;
z := local z;

S := QQ[x_1..x_n];
T := QQ[z_1..z_k];

-- In the ring S, first determine all monomials of degree D. Second,
-- determine which pair of monomials have a LCM of degree D+1

monoS := (monomialIdeal(gens S))^d;
Sd := flatten entries gens monoS;
Pairs := {};
for i from 1 to k do (
for j from (i + 1) to k do (
m := lcm(Sd#(i - 1), Sd#(j - 1));
degm := first degree m;
if degm == d + 1 then Pairs = append(Pairs,{z_i, z_j});
);
);

G := graph(T,Pairs);
L := vertices G;
W := {};

while(#L > 0) do (
v := first L;
if #W == 0 then v = L#(random(0, #L - 1))
else (
-- Select a vertex of minimum degree
dV := apply(L, i -> degreeVertex(G, i));
v = L#(position(dV, i -> i == min dV));
);

-- Add that vertex to the set
W = append(W,v);

-- Remove v and its neighbours from L
L = toList(set(L) - ({v} | neighbors(G,v)));
);

-- For each vertex in the maximal set, find its neighbours and add them to that vertex
newGens := apply(W, v -> (v + sum neighbors(G,v)));
edgeIdealGens := apply(edges G, product);

-- Create a new ideal with these new generators
modifiedIdeal := ideal(edgeIdealGens | newGens);
modifiedIdealDim := time dim modifiedIdeal;
upperBound := modifiedIdealDim + #newGens;

-- If the upper bound matches the cardinality of the maximal set, then it must be
-- the spreading number.
if #W == upperBound then print "Lower bound and upper bound equal.";
return upperBound;
);


-----------------------------------------------------------------------------
-----------------------------------------------------------------------------

- Macaulay2 function to compute an upper bound of the spreading number
--
-- Ben Babcock [bababcoc@lakeheadu.ca]
-------------------------------------------------------------------------------

needsPackage "EdgeIdeals";

isClique = method();
isClique(Graph, List) := (G,c) -> (
for i when i < #c do (
for j from i + 1 to (#c - 1) do (
if not isEdge(G, {c#i, c#j}) then return false;
);
);
return true;
);

-----------------------------------------------------------------------------
-----------------------------------------------------------------------------

spreadingNumberBound = (n,d) -> (
if n == 1 or d == 1 then return 1;
if d == 2 then return n;

k := binomial(d + n - 1, n - 1);
S := QQ[x_1..x_n];
T := QQ[z_1..z_k];

-- In the ring S, first determine all monomials of degree d. Second,
-- determine which pair of monomials have a LCM of degree d + 1
Sd := flatten entries gens((monomialIdeal(gens S))^d);
Pairs := {};
legend := {};
for i from 1 to k do (
legend = append(legend, {z_i, Sd#(i - 1)});

for j from (i + 1) to k do (
m := lcm(Sd#(i - 1), Sd#(j-1));
degm := first degree m;
if degm == d + 1 then Pairs = append(Pairs,{z_i,z_j});
);
);

-- This table will match up generators of T with monomials of S_n(d)
legend = hashTable legend;
inverseLegend := applyPairs(legend, (p,q) -> (q,p));

-- Construct the graph S_n(d)
G := graph(T, Pairs);

-- Partition the vertices of G by the transposition (12...n-1)
remainingVertices := vertices G;
parts := while #remainingVertices > 0 list (
currentVertex := first remainingVertices;
vClass := {currentVertex};
while(true) do (
f := factor legend#(currentVertex);
m := 1;

for k when k < #f do (
i := index f#k#0;
if i < (n - 1) then i = (i + 1) % (n - 1);
m = m*(x_(i + 1)^(f#k#1));
);
m = inverseLegend#m;
currentVertex = m;
if m == first remainingVertices then (
if #vClass == 1 then vClass = append(vClass, m);
break;
);
vClass = append(vClass, m);
);
remainingVertices = toList(set(remainingVertices) - set(vClass));
vClass
);

numPartsCliques := 0;
parts = for i when i < #parts list (
if #(unique parts#i) == 1 then continue
else (
if isClique(G, parts#i) then numPartsCliques = numPartsCliques + 1;
sum parts#i
)
);
edgeIdealGens := apply(edges G, product);

modifiedIdeal := ideal(parts | edgeIdealGens);
modifiedIdealDim := time dim(T/modifiedIdeal);
--print modifiedIdealDim;
return modifiedIdealDim + #parts - numPartsCliques;
);

-- Macaulay2 function for computing an upper bound of the covering number
--
-- v1.1
-- Ben Babcock 
------------------------------------------------------------------------------

needsPackage "EdgeIdeals";

getMinimalCover = (n,d) -> (
if n == 1 or d == 1 then return 1;
if d == 2 then return n;

k := binomial(d + n - 1, n - 1);
S := QQ[x_1..x_n];
T := QQ[z_1..z_k];

-- In the ring S, first determine all monomials of degree d. Second,
-- determine which pair of monomials have a LCM of degree d + 1
Sd := flatten entries gens((monomialIdeal(gens S))^d);
Pairs := {};
legend := {};
for i from 1 to k do (
legend = append(legend, {z_i, Sd#(i - 1)});
for j from (i + 1) to k do (
m := lcm(Sd#(i - 1), Sd#(j-1));
degm := first degree m;
if degm == d + 1 then Pairs = append(Pairs,{z_i,z_j});
);
);


-- This table will match up generators of T with monomials of S_n(d)
legend = hashTable legend;
inverseLegend := applyPairs(legend, (p,q) -> (q,p));

-- Construct the graph S_n(d)
G := graph(T, Pairs);

-- Find all upward cliques (each is identified by a monomial from S_n(d - 1))
SdMinus1 := flatten entries gens((monomialIdeal(gens S))^(d - 1));
upwardCliques := hashTable apply(SdMinus1, i -> ({i, for j from 1 to n list i*x_j}));

-- Our initial set, for the cover generator, will be the cliques containing
-- x_1^d,...,x_n^d, because these each belong to one and only one clique.
upCover := for i from 1 to n list upwardCliques#(x_i^(d - 1));

-- Get a list of all vertices in the cover
W := flatten upCover;
coveredVertices := unique apply(W, w -> inverseLegend#w);

while true do (
uncoveredVertices := toList(set(vertices G) - coveredVertices);
if #uncoveredVertices == 0 then break;

-- Select a random vertex and find the upward cliques associated with it
v := legend#(uncoveredVertices#(random(0, #uncoveredVertices - 1)));
U := select(values upwardCliques, i -> member(v, i));

candidateClique := first U;
for i from 1 when i < #U do (
-- Determine how many vertices are in the cover
vInCover := select(U#i, j -> member(j, W));
if #vInCover == 0 or #vInCover < #select(candidateClique, j -> member(j, W)) then candidateClique = U#i;
);

-- Add the candidate clique to the cover
upCover = append(upCover, candidateClique);
coveredVertices = coveredVertices | apply(candidateClique, w -> inverseLegend#w);
);

-- Minimize the cover
while true do (
-- Get the frequency for each vertex (how many cliques contain it)
vertexFreq := new MutableHashTable;
for i when i < #upCover do (
for j when j < #(upCover#i) do (
v := (upCover#i)#j;
if vertexFreq#?v then vertexFreq#v = vertexFreq#v + 1
else vertexFreq#v = 1;
);
);

-- Iterate over the cliques. If there is a clique in which no vertex is of frequency
-- 1, then the clique is not essential to the cover and may be discarded
discard := for i when i < #upCover do (
c := upCover#i;
essential := false;
for j when j < #c do (
if vertexFreq#(c#j) == 1 then (essential = true; break;);
);
if not essential then break i;
);

if discard =!= null then upCover = drop(upCover, {discard,discard})
else break;
);

return #upCover;
);

--------------------------------------------------------------------------
--------------------------------------------------------------------------

needsPackage "EdgeIdeals";

coveringNumberBound = (n,d) -> (
<< "Degree " << d << " in " << n << " variables " << endl;
if n == 1 or d == 1 then return 1;
if d == 2 then return n;

k := binomial(d + n - 1, n - 1);
S := QQ[x_1..x_n];
T := QQ[z_1..z_k];

sGroup := symmetricGroup S;

-- In the ring S, first determine all monomials of degree d. Second,
-- determine which pair of monomials have a LCM of degree d + 1
Sd := flatten entries gens((monomialIdeal(gens S))^d);
Pairs := {};
legend := {};
for i from 1 to k do (
legend = append(legend, {z_i, Sd#(i - 1)});
for j from (i + 1) to k do (
m := lcm(Sd#(i - 1), Sd#(j-1));
degm := first degree m;
if degm == d + 1 then Pairs = append(Pairs,{z_i,z_j});
);
);


-- This table will match up generators of T with monomials of S_n(d)
legend = hashTable legend;
inverseLegend := applyPairs(legend, (p,q) -> (q,p));

-- Construct the graph S_n(d)
G := graph(T, Pairs);

-- Find all upward cliques (each is identified by a monomial from S_n(d - 1))
SdMinus1 := flatten entries gens((monomialIdeal(gens S))^(d - 1));
upwardCliques := hashTable apply(SdMinus1, i -> ({i, for j from 1 to n list i*x_j}));

-- Our initial set, for the cover generator, will be the cliques containing
-- x_1^d,...,x_n^d, because these each belong to one and only one clique.
upCover := for i from 1 to n list upwardCliques#(x_i^(d - 1));

-- Get a list of all vertices in the cover
W := unique flatten upCover;

graphOrbits := reverse orbits(n,d);
for orb when orb < #graphOrbits do (
--<< "Orbit: " << graphOrbits#orb << endl;
f := factor product for i when i < #(graphOrbits#orb) list x_(i+1)^(graphOrbits#orb#i);
for i when i < #sGroup do (
sigma := sGroup#i;
v := product for j when j < #f list (x_(sigma#(index(f#j#0)) + 1)^(f#j#1));

-- If v is in the cover, skip it
if member(v, W) then continue;

-- Find the upward cliques that contain v
U := select(values upwardCliques, i -> member(v, i));

candidateClique := first U;
for i from 1 when i < #U do (
-- Determine how many vertices are in the cover
vInCover := select(U#i, j -> member(j, W));
if #vInCover == 0 or #vInCover < #select(candidateClique, j -> member(j, W)) then candidateClique = U#i;
);

-- Add the candidate clique to the cover
upCover = append(upCover, candidateClique);
W = W | candidateClique;
);
);

print ("Unminimized cover (size: " | #upCover | ")");
--print upCover;

-- Minimize the cover
while true do (
-- Get the frequency for each vertex (how many cliques contain it)
vertexFreq := new MutableHashTable;
for i when i < #upCover do (
for j when j < #(upCover#i) do (
v := (upCover#i)#j;
if vertexFreq#?v then vertexFreq#v = vertexFreq#v + 1
else vertexFreq#v = 1;
);
);

-- Iterate over the cliques. If there is a clique in which no vertex is of frequency
-- 1, then the clique is not essential to the cover and may be discarded
discard := for i when i < #upCover do (
c := upCover#i;
essential := false;
for j when j < #c do (
if vertexFreq#(c#j) == 1 then (essential = true; break;);
);
if not essential then break i;
);

if discard =!= null then upCover = drop(upCover, {discard,discard})
else break;
);

print("Minimized cover (size " | #upCover | ")");
--print upCover;

return #upCover;
);


-- List all the orbits of S_n(d)
-- These correspond to the partitions of d of length at most n
orbits = method();
orbits(ZZ,ZZ) := (n,d) -> (
-- We want partitions of d
parts := partitions d;

-- Now keep only those partitions of length n or smaller
parts = select(parts, i -> (#i <= n));

-- A plain list is fine, and we want to pad shorter entries
return apply(parts, i -> (
i = toList i;
if #i < n then (
difference := n - #i;
zeros := for j from 1 to difference list 0;
return (i | zeros);
) else return i;
));
);


-- Return a list of elements of the symmetric group of degree n
symmetricGroup = method();
symmetricGroup(Ring) := R -> (
-- This line from the Monomial Ideals chapter, by Serkan Hosten and Gregory Smith,
-- of Computations in Algebraic Geometry with Macaulay2
L := terms det matrix table(numgens R, gens R, (i,r) -> r^(i + 1));
sGroup := {};

-- Construct a hashtable to represent each permutation (starting from 0 instead of 1)
for i from 0 to (#L - 1) do (
elem := {};
f := factor L#i;
for j from 0 to (#f - 1) do (
if index f#j#0 =!= null then elem = append(elem, {index f#j#0, f#j#1 - 1});
);
sGroup = append(sGroup, hashTable elem);
);

return sGroup;
);