Optimizing user defined graphicQ, test of graphicness of a sequence

Question

A sequence of non-negative integers is graphic if it is the degree sequence of some simple graph. Graph realization problem is the decision problem where it is asked whether a given sequence is graphic or not. Graphicness can be tested indirectly using DegreeGraphDistribution, RandomGraph, and SimpleGraphQ, or with the Combinatorica function GraphicQ.

The motivation for a user defined test of graphicness is two-fold: on one hand, with the existing built-ins, there is no direct test afaik, and on the other hand, I don't like loading Combinatorica because of compatibility issues. I have implemented the test of Erdős–Gallai theorem as it is described in Aigner and Triesch 1994 but I think it might still be possible to improve the code.

The input is a list of nonnegative integers, the output is True if the list is graphic, False otherwise. There are 3 main ingredients: conjugate partition, Durfee square, and Erdős-Gallai condition.

Consider a partition $p=(p_1,p_2,...,p_m)$ of integer $m$, ordered non-increasing. Notice that some $p_i$ might be 0. The conjugate partition is $q=(q_1,q_2,...,q_m)$, where again some elements might be 0 and $q_i$ is the number of non-negative elements of $p-i$.

m = 8; (*integer m*)
p = {4, 1, 1, 1, 1, 0, 0, 0}; (*one possible partition*)
q = Table[Count[(p - i), p_ /; p >= 0], {i, m}] (*conjugate partition*)
(*{5, 1, 1, 1, 0, 0, 0, 0}*)

Here is my attempt to make a faster, compiled function for this.

(*the joining at the end is to make do with one less call to Total*)
cp3 = Compile[{{list, _Integer, 1}},
Table[Total[UnitStep[Subtract[list, i]]], {i, 1, First[list]}]~
Join~Table[0, {Length[list] - First[list]}]
];

Durfee square $f(p)$ is the largest "square" in a partition $p$, in other words the largest $i$ for which $p_i\geq i$.

(*one possible way to go about it*)
Last@Pick[Range[m], Thread[p >= Range[m]]]
(*1*)

Here is my attempt to make this as fast as possible.

durfeesquarec = Compile[{{list, _Integer, 1}},
Catch@Do[
 If[list[[i]] < i,
  Throw[Subtract[i, 1]]],
 {i, Length[list]}]
];

Finally, the condition for graphicness can be checked using conjugate partition and Durfee square. Provided that $m$ is even, $p$ is graphic iff $\sum_{i=1}^k(q_i-p_i)\geq k$ for $k=1,...,f(p)$.

(*does not check non-negativity*)
graphicQ[list_?VectorQ] := Module[{p = list},
  (*sum even?*)
  If[OddQ[Total[list]], Return[False]];
  (*ordered non-increasing?*)
  If[! OrderedQ[Reverse[p]], p = Reverse@Sort@p];
  (*finally, test the graphicness*)
  And @@ Thread[(Accumulate[(cp3[p] - p)[[;; durfeesquarec[p]]]] >= 
      Range[durfeesquarec[p]])]
  ]

graphicQ[_] := False

Test the speed.

(*generate some degree sequences*)
testlist = VertexDegree /@ RandomGraph[{30, 100}, 10000];

(*using preloaded built-ins only*)
SimpleGraphQ@RandomGraph[DegreeGraphDistribution[#]] & /@ testlist // 
Tally // AbsoluteTiming
(*{3.76298, {{True, 10000}}}*)

(*using the GraphicQ from Combinatorica*)
Quiet@Block[{$ContextPath}, Needs["Combinatorica`"]];
Combinatorica`GraphicQ /@ testlist // Tally // AbsoluteTiming
(*{8.08287, {{True, 10000}}}*)

(*graphicQ*)
graphicQ /@ testlist // Tally // AbsoluteTiming
(*{0.542627, {{True, 10000}}}*)

I'm interested in knowing if the code can still be improved, mostly in terms of speed. In particular, I'm wondering if the call to TotalAll during the evaluation of cp3 is a sign that compilation wasn't as good as it could have been. If I wanted it to accept also inputs of the form {2.,1.,1.}, what would be a fast way to check and transform such inputs to lists of integers, e.g., {2,1,1}?

Answer 1

You can also directly compile @Szabolcs code without CCompilerDriver,

ErdosM = Compile[{{input, _Integer, 1}},
  Block[{n, w, b, s, c, k, degrees, res = True},
   If[! EvenQ@Total[input], Return[False]];
   degrees = Reverse@Sort[input];
   n = Length[degrees];
   w = n; b = 0; s = 0; c = 0;
   For[k = 1, k <= n, ++k, b += degrees[[k]];
    c += w - 1;
    While[(w > k && degrees[[w]] <= k), 
     s += degrees[[w]];
     c -= k;
     w--];
    If[(b > c + s), res = False; Break[],
     If[w == k, res = True; Break[]]];
    ];
   res
   ], CompilationTarget -> "C"
  ]

res = ErdosM /@ testlist; // AbsoluteTiming

{0.025278, Null}

Counts[res]

<|True -> 10000|>

Answer 2

The IGraph/M package has a fast graphicality test implemented in C.

<< IGraphM`

?IGGraphicalQ

IGGraphicalQ[degrees] tests if degrees is the degree sequence of any simple undirected graph.

IGGraphicalQ[indegrees, outdegrees] tests if indegrees with outdegrees is the degree sequence of any simple directed graph.

res = IGGraphicalQ /@ testlist; // AbsoluteTiming
(* {0.121092, Null} *)

Counts[res]
(* <|True -> 10000|> *)

I also have a straightforward LibraryLink/C implementation based on this paper by Zoltán Király. This is about as fast as it gets.

res = ErdosGallai /@ testlist; // AbsoluteTiming
(* {0.044923, Null} *)

Counts[res]
(* <|True -> 10000|> *)

---

Here's the code:

Needs["CCompilerDriver`"]

src = "
  #include <WolframLibrary.h>

  DLLEXPORT mint WolframLibrary_getVersion(){
    return WolframLibraryVersion;
  }

  DLLEXPORT int WolframLibrary_initialize( WolframLibraryData libData) {
    return 0;
  }

  DLLEXPORT void WolframLibrary_uninitialize( WolframLibraryData libData) {
    return;
  }

  DLLEXPORT int ErdosGallai(WolframLibraryData libData, mint Argc, MArgument *Args, MArgument Res) {
    MTensor degrees;
    const mint *d, *dims;
    mint n, w, b, s, c;
    mbool res;

    degrees = MArgument_getMTensor(Args[0]);
    d = libData->MTensor_getIntegerData(degrees);
    dims = libData->MTensor_getDimensions(degrees);

    n = dims[0];

    w = n; b = 0; s = 0; c = 0;
    for (mint k=1; k <= n; ++k) {
        b += d[k-1];
        c += w-1;
        while (w > k && d[w-1] <= k) {
            s += d[w-1];
            c -= k;
            w--;
        }
        if (b > c+s) {
            res = False;
            break;
        }
        else if (w == k) {
            res = True;
            break; 
        }   
    }

    MArgument_setBoolean(Res, res);

    return LIBRARY_NO_ERROR;
  }
  ";

CreateLibrary[src, "ErdosGallai"]

eg = LibraryFunctionLoad["ErdosGallai", "ErdosGallai", {{Integer, 1, "Constant"}}, True | False];

ErdosGallai[d : {__Integer}] := EvenQ@Total[d] && eg@Reverse@Sort[d]