The New Combinatorica

Sriram V. Pemmaraju, The University of Iowa
Steven S. Skiena, State University of New York, Stony Brook
May 10, 2001

Introduction

Combinatorica is a Standard Mathematica Add-On written in 1989 by Steven Skiena. It has not been updated since. It has about 230 functions for doing computational discrete mathematics. 
Feedback from users, our own expectations of what Combinatorica should be able to do, various advances in graph theory and combinatorics, faster machines, and better versions of Mathematica are factors that have motivated this rewrite. About 80% of the functions have been rewritten and the package now contains about 350 functions. 
The new Combinatorica provides functions for enumerating, selecting, ranking, and unranking various combinatorial objects such as permutations, combinations, integer partitions, set partitions, Young tableaux, trees, and graphs. It also provides functions to generate various classes of graphs and provides functions for all the standard graph algorithms. The specific ways in which the new Combinatorica improves over the old version are as follows:
	- Improved graph data structure, especially tuned for sparse graphs
	- Functions provided for additional topics such as Set Partitions and Polya Theory
	- Better graphics, with graph drawing significantly improved
	- Many functions substantially speeded up
	- Many old functions now have additonal functionality providing users greater flexibility, 		ease of use, and more error checking
	- New graph instances and graph classes
	- Many miscellaneous new functions
Each of these items is examined below in some detail.

Old Combinatorica vs New Combinatorica

Better Graph Data Structure

A graph is now represented as a triple, the first element is an edge list, the second element is the embedding of the vertices, and the third element is optional graphics information. The main difference is that the adjacency matrix representation has been replaced by an edge list representation. The implications of this change are felt throughout the package - in running time improvements, memory savings, increased functionality, and better graph drawings. The package can now work with graphs that are about 50 times larger than graphs that Old Combinatorica could deal with.

g  =  CompleteGraph[4]

g[[1]]

g[[2]]

g[[0]]

g  =  
      SetGraphOptions[g,  VertexColor -> Red,  
        EdgeColor -> Blue]

ShowGraph[g]

g  =  
      SetGraphOptions[
        g,  {{{1, 2},  VertexStyle -> Disc[Large],  
            VertexColor -> Green}}]

ShowGraph[g,  PlotRange -> Large[0.05]]

Length[g]

g[[1]]

g[[2]]

g[[3]]

g[[4]]

This data structure allows us to store graphics information that pertains to the entire graph or to individual elements such as edges and vertices.

For sparse graphs the savings in memory is dramatic. For dense graphs there is no significant difference in memory usage.

g  =  DiscreteMath`OldCombinatorica`GridGraph[20, 20];

ByteCount[g]

g  =  GridGraph[20,  20];  ByteCount[g]

h  =  
      SetGraphOptions[g,  VertexColor -> Red,  
        EdgeColor -> Blue];  
    ShowGraph[h,  Background -> Yellow]

ByteCount[h]

g  = DiscreteMath`OldCombinatorica`RandomGraph[100,   .5];  
    ByteCount[g]

g  =  RandomGraph[100,   .5];  ByteCount[g]

Additional topics

A number of functions relating to Polya theory and set partitions have been added to the new Combinatorica.

Polya Theory

Groups commonly used in Polya theory can now be generated.

?DihedralGroup

DihedralGroup[4]

The cycle index of groups commonly used in Polya theory can now be computed.

?SymmetricGroupIndex

SymmetricGroupIndex[6,  x]

?DihedralGroupIndex

DihedralGroupIndex[10,  y]

Non-isomorphic instances of various combinatorial objects can now be counted and enumerated. First we show examples involving graphs.

?GraphPolynomial

GraphPolynomial[6, x]

?NumberOfGraphs

Table[NumberOfGraphs[i],  {i,  1,  15}] // TableForm

ShowGraphArray[ListGraphs[6, 3],  Background -> LightBlue]

Now we show examples involving "necklaces." In the first example below, we show all the distinct necklaces with 5 beads, colored red and blue.

?ListNecklaces

ListNecklaces[5,  {r,  r,  b,  b,  y},  Dihedral]

ListNecklaces[5,  {r,  r,  b,  b,  y},  Cyclic]

?NumberOfNecklaces

Table[NumberOfNecklaces[i,  2,  Dihedral],  {i,  10,  20}] // 
      TableForm

?NecklacePolynomial

NecklacePolynomial[10,  {R,  B,  G},  Cyclic]

Set Partitions

?SetPartitions

SetPartitions[4]

?KSetPartitions

KSetPartitions[5,  3]

?RankSetPartition

RankSetPartition[ {{1, 2},  {3, 5},  {4}}]

?UnrankSetPartition

UnrankSetPartition[%,  5]

?RandomSetPartition

RandomSetPartition[10]

?RandomKSetPartition

RandomKSetPartition[10,  4]

?BellB

BellB[300]

The following graph g shows 2-block set partitions of {1, 2, 3, 4, 5}  connected by edges whenever a set partition can be obtained from another by deleting an element from a block and inserting it elsewhere. In the graph below that we show a Hamiltonian cycle in g indicating a "gray code" ordering of the 2-block set partitions.

l52 = 
        Map[StringJoin[Map[ToString, #[[1]]]] &, 
          sp52 = KSetPartitions[5, 2]];
ShowGraph[
      g = SetVertexLabels[
          MakeGraph[
            sp52, (MemberQ[{1, 4}, 
                  Sum[Abs[Position[#1, i][[1, 1]] - Position[#2, 
                            i][[1, 1]]], {i, 5}]]) &,  
            Type -> Undirected], l52], 
      PlotRange -> Large[0.2]]

ShowGraph[Highlight[g, {Partition[HamiltonianCycle[g], 2, 1]}, {Red}], 
      PlotRange -> Large[0.1]]

The graph above has not 1, but more than 6000 Hamiltonian cycles!

Length[HamiltonianCycle[g,  All]]

Our implementation of the Stirling number of the second kind is not recursive. It uses an identity that expresses these Stirling numbers as the signed sum of binomial numbers. This makes our implementation faster than the Mathematica implementation.

StirlingSecond[100,  50]

Table[
      Timing[StirlingS2[100 i, 50 i];]/
        Timing[StirlingSecond[100 i, 50 i];], {i, 2, 7}]

Better Graphics

Old Combinatorica cannot display the colored graph. On the other hard, consider whatNew Combinatorica can do.

c = VertexColoring[ g  =  Wheel[20]]

ShowGraph[
      Highlight[ g,  
        Table[ Flatten[Position[c, i]],  {i,  Max[c]}]]]

Here are some pictures we can easily create in the New Combinatorica.

g  =  GraphProduct[GridGraph[2, 2],  Cycle[5]];
ShowGraph[g,  VertexColor -> LightBlue,  EdgeColor -> Yellow,  
      Background -> Black]

g  =  GridGraph[3, 3];
ShowGraph[g,  VertexStyle -> Disc[Large],  VertexNumber -> Center,  
      VertexNumberColor -> Yellow,  Background -> Yellow]

g  =  
        InduceSubgraph[ GridGraph[10, 10],  
          RandomSubset[ Range[100]]];
m  =  BipartiteMatching[ g ]

ShowGraph[Highlight[g,  {m},  {Red}],  
      Background -> LightBlue]

ShowGraph[g  =  Hypercube[5],  VertexColor -> Red]

ShowGraph[g,  VertexColor -> Red,  
      PlotRange -> Zoom[{1, 2, 3, 4}]]

ShowGraphArray[
      Partition[l = Union[Map[(CirculantGraph[9, #]) &, Subsets[8]]], 4, 
        4]]

g  =  Star[7];  ShowGraph[g,  VertexNumber -> On]

ShowGraph[
      AddEdges[g,  {{1,  7},  {1,  1},  {6,  7},  {7,  
            6}}]]

Improvement in Run Times

Many functions in the Old Combinatorica are quite slow. This is for a variety of reasons. Some of the reasons cannot be remedied. For example, the overhead of using Mathematica or the inherent intractability of certain problems. However, using better algorithms, better data structures, better programming practices, and new Mathematica features we have speeded up most functions, some by several orders of magnitude.

The speedup you see in the following examples is because the adjacency matrix representation for graphs has been replaced by an edge list representation. This is ideal for the sparse graphs being generated below.

{Timing[ DiscreteMath`OldCombinatorica`Path[300];],  
      Timing[Path[300];]}

{Timing[ DiscreteMath`OldCombinatorica`RandomTree[200];], 
      Timing[RandomTree[200];]}

{Timing[ DiscreteMath`OldCombinatorica`GridGraph[30, 30];], 
      Timing[GridGraph[30,  30];]}

Various utility functions that convert between different representations of graphs have also been speeded up and this has led to speedup of almost all the graph algorithms in the package.

Observe the remarkable speedup in the function ToAdjacencyLists, shown below. ToAdjacencyLists is used in most graph algorithms. Below we attempt to show the difference in the running times of MinimumSpanningTree in Old Combinatorica and New Combinatorica. The attempt is unsuccessful because MinimumSpanningTree in Old Combinatorica takes more time than we could spend waiting! Some of this difference in running time is due to the difference in the running times of ToAdjacencyLists.

g  =  DiscreteMath`OldCombinatorica`GridGraph[30,  30];  
    Timing[ DiscreteMath`OldCombinatorica`ToAdjacencyLists[
          g];]

g  =  GridGraph[30,  30];  
    Timing[ ToAdjacencyLists[g];]

g  =  SetEdgeWeights[GridGraph[30,  30]]; 
    Timing[ MinimumSpanningTree[g];]

g  =   DiscreteMath`OldCombinatorica`GridGraph[20,  20]; 
    Timing[DiscreteMath`OldCombinatorica`MinimumSpanningTree[
          g];]

We can now compute with fairly large size sparse graphs. Here we construct a random subgraph with about 5000 vertices and edges of a graph with about 10,000 vertices and 40,000 edges. The subgraph has 668 connected components and it took about 6 seconds to compute these.

g  =  GridGraph[100, 100];  
    h  =  InduceSubgraph[g,  RandomSubset[Range[10000]]];

{V[h],  M[h]}

ShowGraph[h,  VertexStyle -> Disc[0]]

Timing[c = ConnectedComponents[h];]

Length[c]

Some speedup is achieved by an improvement in implementation - more careful attention is paid to ensuring that the actual running time of the implementation is the same as what is promised by theory. The examples below show speedups achieved in functions such as LineGraph and VertexColoring.

g  =  DiscreteMath`OldCombinatorica`RandomGraph[100,   .3];  
    Timing[DiscreteMath`OldCombinatorica`LineGraph[g];]

g  =  RandomGraph[100,   .3];  Timing[LineGraph[g];]

g  =  DiscreteMath`OldCombinatorica`Wheel[200];  
    Timing[DiscreteMath`OldCombinatorica`VertexColoring[ 
          g ];]

g  =  Wheel[200];  Timing[ VertexColoring[g]; ]

Speedup in some cases was achieved by compiling carefully selected functions. In the example below, LexicographicPermutations is reimplemented as an iterative (rather than as a recursive function) and then compiled.

Timing[DiscreteMath`OldCombinatorica`LexicographicPermutations[
          Range[8]];]

Timing[LexicographicPermutations[8];]

In some implementations, we took advantage of the new packed array implementation of lists. Here is a comparison between our new implementation of ToCycles and the old implementation of ToCycles from the Permutations package.

p = DiscreteMath`Permutations`RandomPermutation[1000];  
    Timing[DiscreteMath`Permutations`ToCycles[p];]

p  = RandomPermutation[1000];  Timing[ToCycles[p];]

In fact, ToCyles is so fast that the function SignaturePermutation that computes the sign of a permutation is significantly faster than the corresponding Mathematica function Signature.

p  =  RandomPermutation[5000];  {Timing[Signature[p];],  
      Timing[SignaturePermutation[p];]}

Better Functionality

A variety of old functions are now easier to use. A variety of new user friendly functions have been added, Options have been used to make several of the functions more user friendly and flexible. For example, AddEdges allows many variations of edge specifications.

?AddEdges

g  =  Wheel[10];  ShowGraph[g,  VertexNumber -> On]

ShowGraph[
      g = AddEdges[
          g,  {{{2,  10},  EdgeColor -> Red,  
              EdgeStyle -> Fat}}],  VertexNumber -> On]

ShowGraph[g  =  AddEdges[g,  {3,  1}],  
      VertexNumber -> On]

Here is another example, in which we show how the function BreadthFirstTraversal has become much more flexible. First compare the usage messages for the new BreadthFirstTraversal and the old BreadthFirstTraversal. Some additional tags allow the user to get different kinds of information from BreadthFirstTraversal.

?BreadthFirstTraversal

?DiscreteMath`OldCombinatorica`BreadthFirstTraversal

g  =  GridGraph[5, 5];

BreadthFirstTraversal[g,  1]

e = BreadthFirstTraversal[g,  1,  Edge]

ShowGraph[Highlight[g,  {e},  {HotPink}]]

ShowGraphArray[{g,  BreadthFirstTraversal[g, 1, Tree]}]

BreadthFirstTraversal[g,  1,  Level]

Various functions such as Eccentricity, Diameter, TwoColoring, BipartiteQ, and others now  simply call BreadthFirstTraversal with an appropriate argument.

Functions that generate graphs often come with options that allow the user to choose between directed or undirected graphs. Here is an example.

?CompleteGraph

ShowGraph[ CompleteGraph[5,  Type -> Directed] ]

SetOptions[ CompleteGraph,  Type -> Directed]

ShowGraph[ CompleteGraph[5] ]

Functions that translate between different graph representations now have options that allow the user to pay attention to multiple edges and self-loops. In the following example,  the function ToAdjacencyLists, in its default version, makes sure that self-loops and multiple edges show up in the adjacency lists. Using the option Type->Simple, we can force ToAdjacencyLists to ignore the self-loops and multiple edges.

g  =  AddEdges[Star[5],  {{1,  5},  {2,  2}}]

ShowGraph[ g ]

ToAdjacencyLists[g]

ToAdjacencyLists[g,  Type -> Simple]

New Graph Instances and Classes

?GridGraph

ShowGraph[GridGraph[3,  4,  5] ]

?ShuffleExchangeGraph

ShowGraph[ ShuffleExchangeGraph[3] ]

?ButterflyGraph

ShowGraph[ ButterflyGraph[3] ]

?InversionPoset

ShowGraph[ InversionPoset[4] ]

?BooleanAlgebra

ShowGraph[ BooleanAlgebra[4] ]

?CageGraph

ShowGraphArray[{Table[CageGraph[r],  {r,  3,  5}],  
        Table[CageGraph[r],  {r,  6,  8}]}]

Miscellaneous New Functions

Here we show a new function SetEdgeWeights that provides a flexible way of assigning  weights to edges in a graph.

?SetEdgeWeights

g  =  
      SetEdgeWeights[Wheel[10],  WeightingFunction -> RandomInteger,  
        WeightRange -> {3,  5}];  GetEdgeWeights[g]

Options[SetEdgeWeights]

?SetOptions

SetOptions[SetEdgeWeights,  WeightingFunction -> RandomInteger,  
      WeightRange -> {0, 1}]

g  =  SetEdgeWeights[Wheel[10]]; GetEdgeWeights[g]

Here is another example of a useful new function. Shortest paths can now be computed even in the presence of negative edge weights.

?BellmanFord

g  = 
      AddEdges[ 
        Cycle[10,  
          Type -> Directed],  {{4,  1},  {6,  1},  {9,  
            3}}]

ShowGraph[g,  VertexNumber -> On]

g  =  
      SetEdgeWeights[g,  WeightingFunction  ->  RandomInteger,  
        WeightRange -> {-4,  5}]

Edges[g,  EdgeWeight]

BellmanFord[g,  1]

Dijkstra[g,  1]

Another example of a useful new function is TreeIsomorphismQ.

?TreeIsomorphismQ

g  =  RandomTree[40];  
    h  =  PermuteSubgraph[g,  RandomPermutation[40]]

TreeIsomorphismQ[g,  h]

IsomorphicQ[g,  h]

{Timing[TreeIsomorphismQ[g,  h];],  
      Timing[IsomorphicQ[g,  h];]}