1
$\begingroup$

I have some code to generate a random geometric graph in a domain of specificed shape (here an H-shape). Edges cannot exist when they intersect the domain boundary (so we need line of sight for edges to exist).

Here is the domain shape:

    reg[a_, b_, w_] := 
 RegionUnion[Rectangle[{0, 0}, {a, b}], 
  Rectangle[{a, (b - 1)/2}, {a + w, (b + 1)/2}], 
  Rectangle[{a + w, 0}, {a + w + a, b}]](*Define H-shaped region*)

    regcomp[a_, b_, w_] := 
 RegionUnion[Rectangle[{a, 0}, {a + w, (b - 1)/2}], 
  Rectangle[{a, (b + 1)/2}, {a + w, 
    b}]](*Definte its complement*)

    Region@reg[1, 2, 4]
    Region@regcomp[1, 2, 4]

enter image description here

This code uses RegionIntersection to check each edge, and remove it accordingly.

n = 100;
range = 1;
reg[a_, b_, w_] := 
 RegionUnion[Rectangle[{0, 0}, {a, b}], 
  Rectangle[{a, (b - 1)/2}, {a + w, (b + 1)/2}], 
  Rectangle[{a + w, 0}, {a + w + a, b}]](*Define H-shaped region*)

regcomp[a_, b_, w_] := 
 RegionUnion[Rectangle[{a, 0}, {a + w, (b - 1)/2}], 
  Rectangle[{a, (b + 1)/2}, {a + w, 
    b}]](*Definte its complement*)
edgesold = 
 Function[{subsets, b}, 
  Pick[subsets, 
   UnitStep[
    RandomReal[{0, 1}, Length@subsets] - 
     Exp[-b Dot[(subsets[[All, 1]] - subsets[[All, 2]])^2, {1., 
         1.}]]], 0]];(*Make edges for the graph*)
edges[v_] := 
 Module[{e1}, e1 = edgesold[Subsets[v, {2}], 1/range^2]; 
  e2 = Select[e1, 
    RegionDimension@
       RegionIntersection[Line[#], 
        regcomp[1, 4, 2]] > -\[Infinity] &]; 
  DeleteCases[e1, Alternatives @@ e2]
  ] (*Remove edges which intersect boundary*)

gr6[vert_] := 
 Graph[vert, UndirectedEdge @@@ edgesold[Subsets[vert, {2}], beta], 
  VertexCoordinates -> vert];(*Make random geometric graph*)
v = 
 RandomReal[{0, 4}, {n, 2}];
v = Select[v, 
  RegionMember[reg[1, 4, 2], #] == 
    True &];(*Selcted only nodes in H-shaped region*)
graph6 = 
 gr6[v];
rg1 = Graph[v, UndirectedEdge @@@ edges[v], 
  VertexCoordinates -> 
   v](*Build graph on nodes, using only edge not intersecting the \
boundary*)

enter image description here

Is my arrangement of the commands leading to the long time this takes to run? Or, for example, is it better to call RegionIntersection without Select? What's the best way to do this in Mathematica? I need to run this on at least n = 10,000 nodes, which seems to take days to run at the moment.

$\endgroup$

1 Answer 1

3
$\begingroup$

You can speed up by adding more manual checks like convexity and internal intersection code. Here's the example.

regList[a_, b_, w_] := {Rectangle[{0, 0}, {a, b}], 
  Rectangle[{0, (b - 1)/2}, {2 a + w, (b + 1)/2}], 
  Rectangle[{a + w, 0}, {a + w + a, b}]} (* for convexity check *)

hGraph[n_, range_ ,{a_,b_,w_}] :=
    Block[{br, v, epos, pick, memQ, edges, edge1, edge2, member, bounds},
        br = BoundaryDiscretizeRegion[reg[a,b,w], MaxCellMeasure->Infinity];
        v = Select[RandomReal[MinMax[RegionBounds[br]],{n,2}], RegionMember[br,#]==True&];
        epos = Subsets[Range@Length[v], {2}];
        epos = Pick[epos, UnitStep[RandomReal[{0,1}, Length@epos]-Exp[-(1 / range^2) Dot[(v[[epos[[All,1]]]]-v[[epos[[All,2]]]])^2,{1.,1.}]]],0];
        memQ = RegionMember[BoundaryDiscretizeRegion[#]] & /@ regList[1, 4, 2];
        member = Transpose[Through[memQ[v]]];
        pick = AnyTrue[Transpose[member[[#]]], And@@#&]& /@ epos;
        edges = v[[#]] & /@ epos;
        edge1 = Pick[edges, pick];
        edge2 = Pick[edges, pick, False];
        bounds = MeshPrimitives[BoundaryDiscretizeRegion[regcomp[a, b, w], MaxCellMeasure->Infinity],1];
        pick = ParallelMap[Region`Mesh`IntersectionFreeSegmentsQ[{bounds, Line[#]},"Ignore"->{"EndPointsTouching","PointTouchesEndPoint","PointOverlap"}]&, edge2];
        edge2 = Pick[edge2, pick];
        Graph[v,UndirectedEdge@@@Join[edge1, edge2],VertexCoordinates->v]
    ]

In[362]:= g = hGraph[200, 1, {1, 4, 2}]; // AbsoluteTiming

Out[362]= {0.075641, Null}

In[364]:= g = hGraph[10000, 1, {1, 4, 2}]; // AbsoluteTiming

Out[364]= {74.0034, Null}

Show[{Region[Style[reg[1, 4, 2], Opacity[.6]]], 
  GraphPlot[hGraph[600, 1, {1, 4, 2}], VertexShapeFunction -> "Point",
    EdgeShapeFunction -> "Line", EdgeStyle -> Thickness[.0002], 
   BaseStyle -> Black]}]

enter image description here

$\endgroup$
3
  • $\begingroup$ hGraph[600, 1, {1, 4, 2}] seems no output in 13.0.1. $\endgroup$
    – cvgmt
    Apr 19, 2022 at 1:03
  • $\begingroup$ @cvgmt maybe GraphPlot[hGraph[600, 1, {1, 4, 2}]] ? $\endgroup$
    – halmir
    Apr 19, 2022 at 3:18
  • $\begingroup$ Thank you, convexity checks i.e. checking if endpoints are both in a large convex subdomain is what I was trying. Discretising the boundary is also very useful. Also I used ParallelMap when applying the intersection check, which make a big difference on my quad-core Desktop. $\endgroup$
    – apg
    Apr 19, 2022 at 12:44

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.