Optimizing user defined graphicQ, test of graphicness of a sequence

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*)
(*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];

SimpleGraphQ@RandomGraph[DegreeGraphDistribution[#]] & /@ testlist //
Tally // AbsoluteTiming
(*{3.76298, {{True, 10000}}}*)

(*using the GraphicQ from Combinatorica*)
Quiet@Block[{\$ContextPath}, Needs["Combinatorica"]];
CombinatoricaGraphicQ /@ 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}?

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|>

• Thanks for this! As a side note, it's a bit ironic that since I have never really used anything but Mathematica, I'm really bad at writing procedural code which you apparently end up needing when you want the best speed possible.
– Kiro
May 31 '17 at 10:07

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);
d = libData->MTensor_getIntegerData(degrees);
dims = libData->MTensor_getDimensions(degrees);

n = dims;

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]