Stepdis[d_]:=Join[Table[UnitVector[d,k],{k,d}],Table[-UnitVector[d,k],{k,d}]];
T[N_]:=Flatten[Table[{i,j,k},{i,0,N-1},{j,0,N-1},{k,0,N-1}],2]

Step distribution of random walk and 3-dimensional torus of size length N.

TraceRW[x_,p_,N_]:=Mod[Accumulate[Join[{x},RandomChoice[Stepdis[3],p]]],N];

Computes the trace associated to a random walk started at x on the torus T after p step.

palpha[alpha_,N_]:=IntegerPart[alpha*1.516386059151978*Log[N^3]*N^3];

Computes the number of iterations to be done at level alpha (for the alpha late-points as in Miller, Sousi) in a box of size N.

Let us now compute the vacant sets for alpha=0.1, 0.2, ... 1, and save the corresponding plot.

Nu=230;x={0,0,0};
$InitialDirectory
Vu=T[Nu];
RW={};
alpha=0;
Do[alpha=alpha+0.1;
Print[alpha];
RW=TraceRW[x,palpha[0.1,Nu],Nu];
Vu=Complement[Vu,RW];
x=RW[[Nu]];
If[Length[Vu]<2,MindisV=Infinity,ne=Nearest[Vu];
MindisV=Min[MapThread[ManhattanDistance,{Vu,ne[Vu,2][[All,2]]}]]];
Export[StringJoin["Vacant_set_torus_size=",TextString[Nu],"alpha=",TextString[alpha],".jpg"],ListPointPlot3D[Vu,BoxRatios->{1,1,1},AxesStyle->Opacity[0],TicksStyle->Opacity[0],PlotLegends->StringJoin["closest=",TextString[MindisV]]]
Print[]]
,{k,1,10}]

pz[z_,N_]:=IntegerPart[1.516386059151978*(Log[N^3]+z)*N^3];

Computes the number of iterations to be done at level z (as in Belius) in a box of size N .

Let us now compute the vacant sets for z=-1.5, - 0.5, 0.5 and 1.5, and save the corresponding plot.

Nu=180;
x={0,0,0};
$InitialDirectory
z=-2.5;
RW=TraceRW[x,pz[z,Nu],Nu];
Vu=Complement[T[Nu],RW];
Do[z=z+1;
Print[z];
RW=TraceRW[x,pz[z,Nu]-pz[z-1,Nu],Nu];
Vu=Complement[Vu,RW];
x=RW[[Nu]];
If[Length[Vu]<2,MindisV=Infinity,ne=Nearest[Vu];
MindisV=Min[MapThread[ManhattanDistance,{Vu,ne[Vu,2][[All,2]]}]]];
Export[StringJoin["Vacant_set_torus_size=",TextString[Nu],"z=",TextString[z],".jpg"],ListPointPlot3D[Vu,BoxRatios->{1,1,1},AxesStyle->Opacity[0],TicksStyle->Opacity[0],PlotLegends->StringJoin["closest=",TextString[MindisV]]]]
,{k,1,4}]

In order to illustrate the role of alpha_*, let us now compute the vacant set at alpha=0.6, 0.65, 0.7, 0.75 (recall alpha_*=0.67), and save the plot so that the points which have at least 1 neighbor at distance < K=4 are highlighted in red.

Nu=100;x={0,0,0};alpha=0;RW=TraceRW[x,palpha[alpha,Nu],Nu];Vu=Complement[T[Nu],RW];K=4;
Timing[Do[alpha=alpha+0.05;
Print[alpha];
RW=TraceRW[x,palpha[0.05,Nu],Nu];
Vu=Complement[Vu,RW];
x=RW[[Nu]];
If[Length[Vu]<2,MindisV=Infinity,ne=Nearest[Vu];
Distance=MapThread[ManhattanDistance,{Vu,ne[Vu,2][[All,2]]}]];
Neighbor={};
Do[If[Distance[[k]]<=K,Neighbor=Join[Neighbor,{Vu[[k]]}]]
,{k,1,Length[Distance]}];
NoNeighbor=Complement[Vu,Neighbor];
MindisV=Min[Distance];
Export[StringJoin["Vacant_set_torus_size=",TextString[Nu],"alpha=",TextString[alpha],"red_distance_smaller",TextString[K],".jpg"],ListPointPlot3D[{NoNeighbor,Neighbor},BoxRatios->{1,1,1},AxesStyle->Opacity[0],TicksStyle->Opacity[0],PlotLegends->StringJoin["closest=",TextString[MindisV]]],PlotStyle->{Directive[{Blue,PointSize[0.012]}],Directive[{Red,PointSize[0.012]}]}];
Print[ListPointPlot3D[{NoNeighbor,Neighbor},BoxRatios->{1,1,1},AxesStyle->Opacity[0],TicksStyle->Opacity[0],PlotLegends->StringJoin["closest=",TextString[MindisV]],PlotStyle->{Directive[{Blue,PointSize[0.012]}],Directive[{Red,PointSize[0.012]}]}]]
,{k,1,12}]]

Bernoulli[N_,p_]:=RandomInteger[{1,N},{RandomVariate[BinomialDistribution[N^3,p]],3}]

Nu=400;alpha=0.6;
Balpha=Bernoulli[Nu,Nu^(-3*alpha)];

K=1;If[Length[Balpha]<2,MindisV=Infinity,ne=Nearest[Balpha];
Distance=MapThread[ManhattanDistance,{Balpha,ne[Balpha,2][[All,2]]}]];
Neighbor={};
Do[If[Distance[[k]]<=K,Neighbor=Join[Neighbor,{Balpha[[k]]}]]
,{k,1,Length[Distance]}];
NoNeighbor=Complement[Balpha,Neighbor];
MindisV=Min[Distance];
Export[StringJoin["Bernoulli_torus_size=",TextString[Nu],"alpha=",TextString[alpha],"red_distance_smaller",TextString[K],".jpg"],ListPointPlot3D[{NoNeighbor,Neighbor},ImageSize->1000,BoxRatios->{1,1,1},AxesStyle->Opacity[0],TicksStyle->Opacity[0]],PlotStyle->{Directive[{Blue,PointSize[0.012]}],Directive[{Red,PointSize[0.012]}]}];
Print[ListPointPlot3D[{NoNeighbor,Neighbor},ImageSize->1000,BoxRatios->{1,1,1},AxesStyle->Opacity[0],TicksStyle->Opacity[0],PlotLegends->StringJoin["closest=",TextString[MindisV]],PlotStyle->{Directive[{Darker[Green],PointSize[0.008]}],Directive[{Red,PointSize[0.008]}]}]]

Nu=400;x={0,0,0};step=0.05;goal=0.6;alpha=0;RW={};Vu=T[Nu];K=1;
Do[alpha=alpha+step;
RW=TraceRW[x,palpha[step,Nu],Nu];
Vu=Complement[Vu,RW];
x=RW[[Nu]];
,{k,1,goal/step}]

If[Length[Vu]<2,MindisV=Infinity,ne=Nearest[Vu];
Distance=MapThread[ManhattanDistance,{Vu,ne[Vu,2][[All,2]]}]];
Neighbor={};
Do[If[Distance[[k]]<=K,Neighbor=Join[Neighbor,{Vu[[k]]}]]
,{k,1,Length[Distance]}];
NoNeighbor=Complement[Vu,Neighbor];
MindisV=Min[Distance];
Export[StringJoin["Vacant_set_torus_size=",TextString[Nu],"alpha=",TextString[alpha],"red_distance_smaller",TextString[K],".jpg"],ListPointPlot3D[{NoNeighbor,Neighbor},ImageSize->1000,BoxRatios->{1,1,1},AxesStyle->Opacity[0],TicksStyle->Opacity[0],PlotLegends->StringJoin["closest=",TextString[MindisV]]],PlotStyle->{Directive[{Blue,PointSize[0.01]}],Directive[{Red,PointSize[0.01]}]}];
Print[ListPointPlot3D[{NoNeighbor,Neighbor},ImageSize->1000,BoxRatios->{1,1,1},AxesStyle->Opacity[0],TicksStyle->Opacity[0],PlotLegends->StringJoin["closest=",TextString[MindisV]],PlotStyle->{Directive[{Darker[Green],PointSize[0.008]}],Directive[{Red,PointSize[0.008]}]}]]

Nu=650;x={0,0,0};step=0.05;alpha=0;RW={};Vu=T[Nu];K=1;
Do[alpha=alpha+step;Print[alpha];
RW=TraceRW[x,palpha[step,Nu],Nu];
Vu=Complement[Vu,RW];
x=RW[[Nu]];
,{k,1,12}]
step=0.005;
Do[alpha=alpha+step;Print[alpha];
RW=TraceRW[x,palpha[step,Nu],Nu];
Vu=Complement[Vu,RW];
x=RW[[Nu]];
If[Length[Vu]<2,MindisV=Infinity,ne=Nearest[Vu];
Distance=MapThread[ManhattanDistance,{Vu,ne[Vu,2][[All,2]]}]];
Neighbor={};
Do[If[Distance[[k]]<=K,Neighbor=Join[Neighbor,{Vu[[k]]}]]
,{k,1,Length[Distance]}];
NoNeighbor=Complement[Vu,Neighbor];
MindisV=Min[Distance];
Print[ListPointPlot3D[{NoNeighbor,Neighbor},ImageSize->1000,BoxRatios->{1,1,1},AxesStyle->Opacity[0],TicksStyle->Opacity[0],PlotLegends->StringJoin["closest=",TextString[MindisV]],PlotStyle->{Directive[{Blue,PointSize[0.007]}],Directive[{Red,PointSize[0.007]}]}]]
,{k,1,30}]