Random Geometric Graphs, but fast

Question

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?

Answer 1

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 *)
     Developer`ToPackedArray[{{-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.

Answer 2

CompiledFunctionTools`CompilePrint@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 = Compile`GetElement[vert, i];
  Table[
   {Compile`GetElement[vert, j], a},
   {j, 1, i-1}],
  {i, 2, Lenght[vert]}]
 ]

Answer 3

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