Sunday, November 25, 2018

Visualizing paths in configuration space

The goal of this post is to visualize how point configurations induce persistent homology, and how paths between point samples induce changes in the simplicial complexes producing the homology. We use the Čech simplicial complex construction of a finite subset of $\R^N$.

Definition: For $M$ a Riemannian mandifold, $\Conf_n(M):= \{P\subseteq M : |P|=n\}$ is the configuration space of $n$ points on $M$.

The space $\Conf_n(M)$ is itself a topological space, with topology induced by the Hausdorff distance of subsets. Let $\SC$ be the set of abstract simplicial complexes $(V,S)$, where $V$ is a set and $S\subseteq P(V)$ closed under subsets. Let $\uSC$ be the set of unlabeled abstract simplicial complexes, with the natural projection map $\SC\to \uSC$.

Definition: The Čech map is the function $\check C\colon \Conf_n(M)\times \R_{\geqslant 0}\to \SC$ given by $V(\check C(P,r))=P$ and $P'\in S(\check C(P,r))$ whenever $\bigcap_{p\in P'} B(p,r) \neq \emptyset$, for every $P'\subseteq P$. The unlabeled Čech map is the composition of $\check C$ with the projection to $\uSC$.

We will consider the case $M=\R^2$ and $n=4$. To describe an implementation of the Čech map, we only need to consider double and triple intersections. Finding if $B(P_1,r)\cap B(P_2,r)$ is empty or not is easy, but to determine if $B(P_1,r)\cap B(P_2,r)\cap B(P_3,r)$ is empty or not requires more care. Below is an implementation in Mathematica.

(* CechPt : Finds the coordinate where balls of the same radii around a,b,c will first intersect *)
(* Input : 3 coordinates {x, y}. Output : 1 coordinate {x, y} *)
CechPt[a_,b_,c_] := Module[{
    cenx = Det[{{Norm[a]^2, a[[2]], 1}, {Norm[b]^2, b[[2]], 1}, {Norm[c]^2, c[[2]], 1}}],
    ceny = Det[{{a[[1]], Norm[a]^2, 1}, {b[[1]], Norm[b]^2, 1}, {c[[1]], Norm[c]^2, 1}}],
    scal = 2*Det[{{a[[1]], a[[2]], 1}, {b[[1]], b[[2]], 1}, {c[[1]], c[[2]], 1}}]},
   cen = {cenx/scal, ceny/scal};
   If[Max[ArcCos[(b-a).(c-a)/(Norm[b-a]*Norm[c-a])],
           ArcCos[(a-b).(c-b)/(Norm[a-b]*Norm[c-b])],
           ArcCos[(a-c).(b-c)/(Norm[a-c]*Norm[b-c])]] < Pi/2, cen,
     If[Norm[cen-(a+b)/2] < Norm[cen-(a+c)/2],
       If[Norm[cen-(a+b)/2] < Norm[cen-(b+c)/2], (a+b)/2, (b+c)/2],
       If[Norm[cen-(a+c)/2] < Norm[cen-(b+c)/2], (a+c)/2, (b+c)/2]]]];

Here cen is the circumcenter of the input points, which corresponds to our desired point only if it lies within the convex hull of the points. Now $B(P_1,r)\cap B(P_2,r)\cap B(P_3,r)$ is non-empty if and only if the distance from each of $P_1$, $P_2$, $P_3$ to CechPt[$P_1$, $P_2$, $P_3$] is less than or equal to $r$.

Let $\gamma\colon I\to \Conf_4(\R^2)$ be a path, and $\gamma(0)= \{P_1,P_2,P_3,P_4\}$. At each $t\in I$ and for every pair and triple $P'\subseteq\gamma(t)$, we can find the smallest $r$ such that $\bigcap_{p\in P'}B(P,r)\neq \emptyset$. This gives 6 curves for the pairs $P'$, and 4 curves for the triples $P'$, which we can plot all together in Mathematica.

PList[t_] := {P1[t],P2[t],P3[t],P4[t]};
(* Graphs of pairwise distances *)
DistGraph1 = Plot[Table[Norm[pair[[1]]-pair[[2]]]/2, {pair,Subsets[PList[t],{2}]}], 
               {t, 0, 1}, PlotRange -> {{0,1},{0,1.5}}, PlotStyle -> {Gray}, AspectRatio -> 1];
(* Graphs of minimum distance from every triple to its CechPt*)
DistGraph2 = Plot[Table[Max[
                  Table[Norm[triple[[k]]-CechPt@@triple],{k,1,3}]], {triple,Subsets[PList[t],{3}]}],
               {t, 0, 1}, PlotRange -> {{0,1},{0,1.5}}, PlotStyle -> {Orange}, AspectRatio -> 1];
The code is given so that it may be easily generalized to more than 4 points. Next, use the Manipulate command to add interactivity to the graphs.

Manipulate[{
  Show[DistGraph2, DistGraph1],
  Show[
    ParametricPlot[PList[t],{t,0,X[[1]]},PlotRange -> {{-2,2},{-2,2}},PlotStyle -> {Black}],
    Graphics[Join[
      {Opacity[.2],Red}, Table[Disk[point,X[[2]]],{point,PList[X[[1]]]}],
      {Opacity[1],Red}, Table[Circle[point,X[[2]]],{point,PList[X[[1]]]}],
      {Red,Disk[P1[X[[1]]],.05]},
      {Blue,Disk[P2[X[[1]]],.05]},
      {Darker[Green],Disk[P3[X[[1]]],.05]},
      {Yellow,Disk[P4[X[[1]]],.05]}]]],
  Graphics[Join[
    {Black, Thick},
    Flatten[Table[{{Opacity[0], Opacity[.3]}[[Boole[X[[2]] >= Norm[pair[[1]][[1]][X[[1]]] 
               - pair[[2]][[1]][X[[1]]]]/2] + 1]], Line[#[[2]]&/@pair]},
               {pair,Subsets[{{P1,{0,0}},{P2,{2,0}},{P3,{0,2}},{P4,{2,2}}},{2}]}]],
    Flatten[Table[{{Opacity[0], Opacity[.3]}[[Boole[X[[2]] >= Max[Table[Norm[triple[[k]][[1]][X[[1]]]
               - CechPt@@(#[[1]][X[[1]]]&/@triple)], {k,1,3}]]] + 1]], Polygon[#[[2]]&/@triple]},
               {triple,Subsets[{{P1,{0,0}},{P2,{2,0}},{P3,{0,2}},{P4,{2,2}}},{3}]}]],
    {Opacity[1], Red, Disk[{0,0},.07], Blue, Disk[{2,0},.07], Darker[Green], Disk[{0,2},.07],
               Yellow, Disk[{2,2},.07]}]]
  }, {{X, {.1, .1}}, Locator}]
This produces the interactive visualization below, allowing the user to drag the crosshairs on the graph on the left (graphs of when double and triple intersections are reached). The paths of the individual points $P_1,P_2,P_3,P_4$ are in the middle and the image of the unlabeled Čech map is on the right.

The graphs on the left stratify the strip $I\times \R_{\geqslant 0}$, so that the unlabeled Čech map is constant on each stratum. Computing the Betti numbers of each simplicial complex gives the CROCKER plot (see TZH) of the stratified space. We use the Čech instead of the Rips complex, so perhaps this should be called the CROCKEČ plot. The stratified space, 0-dimensional, and 1-dimensional plots are given below.

Here the Betti numbers were computed by inspection, since the complexes are so small. An extension would be to make this computation automatic once the input path $\gamma$ is given.

The Mathematica code for this post is available online.

References: Topaz, Ziegelmeier, Halverson (Topological Data Analysis of Biological Aggregation Models)