# Random Geometric Graphs, but fast

I am optimising the random geometric graph generation ideas in this post. I need to generate a large number of random geometric graphs with probabilistic connectivity (soft random geometric graphs), then sample the path lengths between two specified vertices in order to find their distribution.

I am trying to compile the edge selection functions, as you can see below:

beta = 1.5;
density = 10;
r = 2;
Ngraphs = 80;
edgescompiled =
Compile[{{vert, _Real}, {b, _Real}},
Select[Subsets[vert, {2}],
RandomReal[{0, 1}] < Exp[-b (Norm[#])^2] &],
CompilationTarget -> "C", RuntimeAttributes -> {Listable},
Parallelization -> True];
gr[vert_] :=
Graph[#[[1]] <-> #[[2]] & /@ edgescompiled[vert, beta]] ;
n = RandomVariate[PoissonDistribution[4 density r^2],
Ngraphs]; // AbsoluteTiming
v = Table[
Join[Table[{RandomReal[{-r, r}], RandomReal[{-r, r}]}, {k, 1,
n[[k]]}], {{-r/2, 0}, {r/2, 0}}], {k, 1,
Ngraphs}]; // AbsoluteTiming
graphs = gr[#] & /@ v; // AbsoluteTiming


But I get an error related to Instruction 3 in CompiledFunction. Does this not work because Select is not compilable? If it isn't, is there some other way I can performance tune the code?

There are two parts to this answer. The first leverages Mathematica's built-in vectorized functions for a order-of-magnitude gain in speed, when parallelized, over the OP's self-answer. The second shows how to fix the OP's original problem with Compile, but since it's slow, it's perhaps of only academic interest.

### The OP's answer with timings on my machine

beta = 1.5;
density = 3;
r = 2;
Ngraphs = 1000;
edgescompiled = Compile[{{subsets, _Real, 3}, {b, _Real}},
Select[subsets, RandomReal[{0, 1}] < Exp[-b (Norm[#])^2] &],
CompilationTarget -> "C", RuntimeAttributes -> {Listable},
Parallelization -> True];
gr[vert_] := Graph[UndirectedEdge @@@ edgescompiled[Subsets[vert, {2}], beta]];
n = RandomVariate[PoissonDistribution[4 density r^2], Ngraphs]; // AbsoluteTiming
v = Table[
Join[{{-r/2, 0}, {r/2, 0}},
Table[{RandomReal[{-r, r}], RandomReal[{-r, r}]}, {k, 1,
n[[k]]}]], {k, 1, Ngraphs}]; // AbsoluteTiming
graphs = gr[#] & /@ v; // AbsoluteTiming
(*
{0.000094, Null}
{0.085466, Null}
{4.68304, Null}
*)


### My solution

We use Pick and machine-precision functions to implement the boolean selection operation (see also this answer, sects. 2 and 2.3). We also construct v as a list of packed arrays.

edges6 = Function[{subsets, b},
Pick[
subsets,
UnitStep[
RandomReal[{0, 1}, Length@subsets] -
Exp[-b Dot[(subsets[[All, 1]] - subsets[[All, 2]])^2, {1., 1.}]]],
0]
];
gr6[vert_] := Graph[vert, UndirectedEdge @@@ edges6[Subsets[vert, {2}], beta],
VertexCoordinates -> vert];

n = RandomVariate[PoissonDistribution[4 density r^2], Ngraphs]; // AbsoluteTiming
(* construct v as a list of packed arrays *)
v = Table[
Join[RandomReal[{-r, r}, {n[[k]], 2}],  (* leverage built-in array generation *)
DeveloperToPackedArray[{{-r/2, 0}, {r/2, 0}}, Real]], {k, 1,
Ngraphs}]; // AbsoluteTiming
graph6 = gr6 /@ v; // AbsoluteTiming
(*
{0.000094, Null}
{0.007307, Null}
{0.687268, Null}
*)


Parallelization helps, but only a little, because machine arithmetic on vector arguments is already parallelized by the MKL. I have 8 virtual cores, but only 4 real cores for doing machine arithmetic. (Frankly, I wouldn't have been surprised at no speed up, but Subsets is probably expensive and not internally parallelized. Perhaps that's why it is a bit faster.)

LaunchKernels[] (* do once only *)

ParallelMap[gr6, v]; // AbsoluteTiming
(*  {0.472053, Null}  *)


### Fixing the OP's Compile

The problem is that Compile does not know the return-type of Subsets. You can specify it with a third argument:

Compile[
{{vert, _Real, 2}, {b, _Real}},
Select[Subsets[vert, {2}], RandomReal[{0, 1}] < Exp[-b (Norm[#])^2] &],
{{Subsets[_, _], _Real, 3}},  (* return type *)
CompilationTarget -> "C"
(*,RuntimeAttributes->{Listable}, Parallelization->True*)
];


This results in a single MainEvaluate (call back to the main kernel) to evaluate Subsets at the beginning. Avoiding MainEvaluate in loops is important; here it is less so. However, the compiled function with MainEvaluate/Subsets cannot be parallelized in the WVM (the environment in which compiled functions are executed). Thus it does not help as much as the other solutions.

• That is a real move forward in my ability to simulate these graphs at large density. Thank you. Also for suggesting how to fix the original compile problem + the rank of the input tensor. – Alexander Kartun-Giles Feb 27 '18 at 14:35

CompiledFunctionToolsCompilePrint@edgescompiled tells me that the problem lies in Subsets not being compilable.

Either you pass the subsets as arguments or you use something like this instead:

Block[{a},
Table[
a = CompileGetElement[vert, i];
Table[
{CompileGetElement[vert, j], a},
{j, 1, i-1}],
{i, 2, Lenght[vert]}]
]

• So maybe I can just pass subsets into the function as a list? – Alexander Kartun-Giles Feb 27 '18 at 12:07
• I will try this also. – Alexander Kartun-Giles Feb 27 '18 at 12:12
• I even think this can be vectorized in a way that the subsets don't have to be built and no Compile is needed. Something like UpperTriangularize[With[{bla = vert^2}, Outer[Plus, bla, bla, 1] ], 1]. I have to leave for the moment. Ping me later if you have further questions. – Henrik Schumacher Feb 27 '18 at 12:15
• Let me test that. Also, I added the missing 1 to {subsets, _Real, 1}, and it started working. – Alexander Kartun-Giles Feb 27 '18 at 12:20

To expand on Hendrik Schumacher's answer, one can pass in the subsets as a list, as such:

beta = 2;
density = 5;
r = 2;
Ngraphs = 100;
edgescompiled =
Compile[{{subsets, _Real, 3}, {b, _Real}},
Select[subsets, RandomReal[{0, 1}] < Exp[-b (Norm[#])^2] &],
CompilationTarget -> "C", RuntimeAttributes -> {Listable},
Parallelization -> True];
gr[vert_] :=
Graph[vert,
UndirectedEdge @@@ edgescompiled[Subsets[vert, {2}], beta]] ;
n = RandomVariate[PoissonDistribution[4 density r^2],
Ngraphs]; // AbsoluteTiming
v = Table[
Join[{{-r, 0}, {r, 0}},
Table[{RandomReal[{-r, r}], RandomReal[{-r, r}]}, {k, 1,
n[[k]]}]], {k, 1, Ngraphs}]; // AbsoluteTiming
graphs = gr[#] & /@ v; // AbsoluteTiming

l = (Length[FindShortestPath[#, {-r, 0}, {r, 0}]] - 1) & /@
graphs; // AbsoluteTiming


and get

{0.000239, Null}

{0.018778, Null}

{0.839959, Null}

{0.000997, Null}

and graphs is built correctly. The path length finding step is many orders of magnitude faster than the graph generation, hence the need for optimisation at this point. Sampling the graph distance between vertex pairs many times for each graph, rather than only once, can also speed this bit up

• You get a slight improvement with gr[vert_] := Graph[UndirectedEdge @@@ edgescompiled[Subsets[vert, {2}], beta]]; – Henrik Schumacher Feb 27 '18 at 12:56
• @HenrikSchumacher yes that is about 25% faster – Alexander Kartun-Giles Feb 27 '18 at 12:58
• With {subsets, _Real, 1}, the variable subsets will be a list (i.e. rank 1) of real numbers, i.e., the coordinates of a single vertex. Is that what it's supposed to be? I thought it was supposed to be a list of pairs of vertices. A list of pairs of vertices has rank 3, not 1. – Michael E2 Feb 27 '18 at 13:19
• I was under the impression it fed in one element at a time. Obviously it is the rank of the input tensor. – Alexander Kartun-Giles Feb 27 '18 at 13:24
• You're welcome. Please see my answer for another approach. – Michael E2 Feb 27 '18 at 14:31