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.