Vietoris-Rips and Lazy Witness Complexes

This Mathematica demonstration, created by Henry Adams and Jan Segert, computes and displays Vietoris-Rips and lazy witness simplicial complexes.

To run this code, you need to click “Enable Dynamics” or select the menu option ”Evaluation -> Dynamic Updating Enabled”. Click and drag on the “filtration parameter” toggle to change the ball size.

WitnessInitialRadii[distanceMatrix_,nu_]:=
Piecewise[{{ConstantArray[0,Dimensions[distanceMatrix][[2]]],nu==0}},Table[Sort[distanceMatrix[[All,i]]][[nu]],{i,Dimensions[distanceMatrix][[2]]}]]

WitnessEdgeMatrix[distanceMatrix_,initialRadii_]:=
Module[{n,m,i,j,edgeMatrix,highlightMatrix},
n=Dimensions[distanceMatrix][[1]];
m=Dimensions[distanceMatrix][[2]];
edgeMatrix = ConstantArray[0,{n,n}];
highlightMatrix = ConstantArray[0,{n,n}];
For[i=1, i≤n, i=i+1,For[j=i+1, j≤n, j=j+1,
edgeMatrix[[i,j]]=Min[Table[Max[{distanceMatrix[[i,k]],distanceMatrix[[j,k]]}- initialRadii[[k]]],{k,m}]];highlightMatrix[[i,j]]=Nearest[Table[Max[{distanceMatrix[[i,k]],distanceMatrix[[j,k]]}- initialRadii[[k]]],{k,m}]->Automatic,0];];];Return[{edgeMatrix,highlightMatrix}];]

DrawBalls[ballCenters_,initialRadii_,halfRadius_,t_]:=Table[{Hue[0.55,1,.95,.6],Disk[ballCenters[[i]],If[halfRadius,(initialRadii[[i]]+t)/2,initialRadii[[i]]+t]]},{i,Length[ballCenters]} ]

DrawHighlights[ballCenters_,initialRadii_,edgeMatrix_,highlightMatrix_,highlight_,halfRadius_,t_]:=If[highlight,Union[Flatten[Table[If[0≤edgeMatrix[[i,j]]-t≤0.1,Union[Table[{Hue[0.6],Circle[ballCenters[[highlightMatrix[[i,j]][[k]]]],If[halfRadius,(initialRadii[[highlightMatrix[[i,j]][[k]]]]+t)/2,initialRadii[[highlightMatrix[[i,j]][[k]]]]+t]]},{k,Length[highlightMatrix[[i,j]]]}]],{}],{i,Length[edgeMatrix]},{j,i+1,Length[edgeMatrix]} ],1]],{}]

DrawData[data_]:=Table[{Hue[0.15,0.9,0.9],Disk[data[[i]], 0.1]},{i,Length[data]} ]

DrawOneSimplices[vertices_,edgeMatrix_,t_]:=Union[Flatten[Table[ If[i<j && edgeMatrix[[i,j]]≤ t,{AbsoluteThickness[1.8], Line[{vertices[[i]],vertices[[j]]}]},Disk[vertices[[i]], 0.1]],{i,Length[vertices]},{j,i,Length[vertices]}]]]

DrawTwoSimplices[vertices_,edgeMatrix_,t_]:=Union[Flatten[Table[ If[edgeMatrix[[i,j]]≤t &&edgeMatrix[[j,k]]≤t &&edgeMatrix[[i,k]]≤t,{Hue[1,.9,.6,.6],EdgeForm[Black],Polygon[{vertices[[i]],vertices[[j]],vertices[[k]]}] },{}],{i,Length[vertices]},{j,i+1,Length[vertices]}, {k,j+1,Length[vertices]}],2]]

Manipulate[
Module[{vertices,n,ballCenters,initialRadii,edgeMatrix,highlightMatrix,nu,distanceMatrix},
Which[
complex=="Rips(data)"∨complex=="Rips(landmarks)",
	Which[complex=="Rips(data)",vertices=data,complex=="Rips(landmarks)",vertices=landmarks];
	n=Length[vertices];
	ballCenters=vertices;
	initialRadii=ConstantArray[0,n];
	edgeMatrix=Table[EuclideanDistance[vertices[[i]],vertices[[j]]],{i,Length[vertices]},{j, Length[vertices]}];
	highlightMatrix=Table[{i,j},{i,n},{j,n}],
complex=="Witness nu=0"∨complex=="Witness nu=1"∨complex=="Witness nu=2",
	Which[complex=="Witness nu=0",nu=0,complex=="Witness nu=1",nu=1,complex=="Witness nu=2",nu=2];
	halfRadius=False;
	vertices=landmarks;
	ballCenters=data;
	distanceMatrix=Table[EuclideanDistance[landmarks[[i]],data[[j]]],{i,Length[landmarks]},{j, Length[data]}];
	initialRadii=WitnessInitialRadii[distanceMatrix, nu];
	{edgeMatrix,highlightMatrix}=WitnessEdgeMatrix[distanceMatrix, initialRadii]];
Graphics[{
	DrawBalls[ballCenters,initialRadii,halfRadius,t],
	DrawHighlights[ballCenters,initialRadii,edgeMatrix,highlightMatrix,highlight,halfRadius,t],
	DrawTwoSimplices[vertices,edgeMatrix, t],
	DrawOneSimplices[vertices,edgeMatrix, t],
	DrawData[Complement[data,vertices]]},
	PlotRange->{{-3,4},{-3,3}},PlotRangeClipping->True,
	ImageSize->{800,500}]
],
Delimiter,Style["Vietoris-Rips and Lazy Witness Complexes",12,Bold],
{{complex,"Rips(data)","Simplicial complex"},{"Rips(data)","Rips(landmarks)","Witness nu=0","Witness nu=1","Witness nu=2"}},
{{highlight,False,"Highlight active balls"},{True,False}},
{{halfRadius,False,"Half radius (Rips only)"},{True,False}},
{{t,0.3,"Filtration parameter "},0,3},
Initialization->(
landmarks = Table[ {2 Cos[theta]+1.4Random[],Sin[theta]+.3Random[]},{theta,0,2 Pi*(6/7),2Pi*(1/7)}];
data = Union[landmarks,Table[ {2 Cos[theta]+1.4Random[],Sin[theta]+.3Random[]},{theta,0,2 Pi*(9/10),2Pi*(1/10)}]];),
SaveDefinitions->True]

To cite this demo:
@software{FiltrationDemos,
  author = {Adams, Henry and Segert, Jan},
  title = {Simplicial complex filtration demonstrations in {M}athematica},
  url = {https://www.math.colostate.edu/~adams/research/},
  date = {2011},
}

XXXX

Contributed by: XXXX