Skip to content

Instantly share code, notes, and snippets.

@bmershon
Last active October 12, 2017 18:48
Show Gist options
  • Star 3 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save bmershon/41bc67cfedf95f7d196d to your computer and use it in GitHub Desktop.
Save bmershon/41bc67cfedf95f7d196d to your computer and use it in GitHub Desktop.
Vietoris-Rips Complex (simplicial complex)

Copyright (c) 2014 Brooks Mershon

The MIT License - http://opensource.org/licenses/MIT

Please excuse the awful coding practices that may be found here.

This was more or less my first run at using D3 and JavaScript in an undergraduate math course (where writing code wasn't even expected, unless it was MATLAB). This was written before modular design, functional programming, and maturity had gripped the author to even a modest extent. This gadget (application seems too strong) was developed as if it were one colossal whiteboard exercise which permitted only small erasures and steady contributions over half of a semester.

The urge to go back and rearrange or modify in any way is strong, but I am choosing to leave this relic of sophomore year as a happy and messy memory of an assignment of which I am quite proud.

I love what I built, but I regret how it was built. I didn't know better.

A final paper was written about the design process of this and other gadgets created for MATH 412, taught by Paul Bendich during the Fall 2014 semester.

Rips-Complex

This gadget demonstrates an application of the Rips-Complex to point-clouds sampled from continuous curves in the plane. The Rips-Complex allows us to calculate the number of connected components as well as the number of 1-cycles (loops) which appear at a given radius.

Points are sampled from an underlying curve with Gaussian noise applied. Persistent homology is computed by reducing the boundary matrix formed by the simplicial complex (at a given radius).

In order to add edges and triangles, a Delaunay triangulation using the existing vertices (points in the point clound) is calculated. Edges are added when the "balls" around the two vertices bounding an edge touch and the edge appears in the Delaunay triangulation.

Triangles are added immediately when their three 1-dimensional faces are present in the simplicial complex. As a result, what you see when the radius is scrubbed all the way to infinity is exactly the Delaunay triangulation of the point cloud.

The purpose of this gadget is to understand how Betti numbers change along the way as we increase the radius from zero to infinity.

The color of a triangle corresponds to the radius at which it appears. Note that the persistence of each 1-cycle formed can be encoded in a 1-dimensional persistence diagram. Such a diagram is replaced in this example by encoding the 1-dimensional persistence with a color scale, as the legend shows.

name d
figure eight "M100,50C0,100,150,350,300,350S300,100,400,50S650,250,500,250S200,0,100,50"
loop "M106,109C68.05,155.05,135.4,356.25,166,375S257.05,252.9,310,234S502.65,273.9,519,249S480.95,89,419,68S143.95,62.949999999999996,106,109"
squares "M260.0028381347656,269.0028381347656V76.00283813476562H102.00283813476562V215.00283813476562H310.0028381347656V143.00283813476562H425.0028381347656V401.0028381347656H318.0028381347656V270.0028381347656H482.0028381347656"
pinched "M107.00283813476562,79.36221313476562C91.85283813476562,74.86221313476562,376.2528381347656,133.51221313476563,374.0028381347656,166.36221313476562S109.85283813476562,265.3622131347656,92.00283813476562,298.3622131347656S197.5528381347656,401.66221313476564,255.00283813476562,386.3622131347656S497.2028381347656,242.41221313476564,475.0028381347656,196.36221313476562S122.15283813476563,83.86221313476562,107.00283813476562,79.36221313476562"
polygon "M165.00283813476562,292.0028381347656L150.00283813476562,84.00283813476562L347.0028381347656,83.00283813476562L562.0028686523438,170.00283813476562L434.0028381347656,364.0028381347656L277.0028381347656,420.0028381347656L68.00283813476562,414.0028381347656L67.00283813476562,205.00283813476562L319.0028381347656,197.00283813476562"
torn "M110.00283813476562,237.72158813476562L129.83617146809894,219.88825480143228C149.66950480143228,202.05492146809894,189.33617146809894,166.38825480143228,213.66950480143228,196.05492146809894C238.00283813476562,225.72158813476562,247.00283813476562,320.7215881347656,263.83617146809894,335.38825480143225C280.66950480143225,350.05492146809894,305.33617146809894,284.38825480143225,314.0028381347656,294.88825480143225C322.66950480143225,305.38825480143225,315.33617146809894,392.05492146809894,365.6695098876953,374.38825480143225C416.00284830729163,356.7215881347656,524.0028584798176,234.72158813476562,546.6695251464844,161.2215887705485C569.336191813151,87.72158940633138,506.6695149739583,62.72159067789713,455.83617655436194,71.55492401123047C405.0028381347656,80.3882573445638,366.0028381347656,123.0549227396647,327.6695048014323,118.88825591405232C289.33617146809894,114.72158908843993,251.66950480143228,63.72159004211425,239.16950480143228,57.388256708780915C226.66950480143228,51.054923375447586,239.33617146809894,89.3882557551066,233.83617146809894,101.55492194493611C228.33617146809894,113.72158813476562,204.66950480143228,99.72158813476562,191.33617146809894,101.22158813476561C178.00283813476562,102.72158813476562,175.00283813476562,119.72158813476562,165.83617146809894,133.38825480143228C156.66950480143228,147.05492146809894,141.33617146809894,157.38825480143228,140.33617146809894,161.8882548014323C139.33617146809894,166.38825480143228,152.66950480143228,165.05492146809894,148.83617146809894,171.55492146809894C145.00283813476562,178.05492146809894,124.00283813476561,192.38825480143228,113.5028381347656,199.55492146809897L103.00283813476562,206.72158813476562"
<!DOCTYPE html>
<meta charset="utf-8">
<style>
.bold_text {
font: bold 48px monospace;
}
.ball {
fill: steelblue;
stroke-width: 0;
stroke: none;
fill-opacity: 0.1;
}
.dot {
fill: black;
stroke-width: 0;
stroke: none;
fill-opacity: 1.0;
}
.edge {
fill: none;
stroke-width: 0.5;
stroke: black;
fill-opacity: 0.3;
}
form {
margin-left: 40px;
margin-top: 20px;
padding-right: 40px;
margin-bottom: 20px;
}
button {
margin-left: 20px;
margin-right:20px;
}
</style>
<body>
<form oninput="output.value = input.value" style="top:10px;left:10px;">
<input type="checkbox" id="curve_checkbox" value="show">Show Sampled Curve
<button id= "sample" type='button'>Sample</button>
<input id="noise-number" type="range" min="0" max="50" value="20" style="width:50;">
<label>Noise</label>
<button id= "betti1_button" type='button'>Compute &#946;<sub>1</sub></button>
<br> <br>
<input id="sample-number" type="range" min="0" max="55" value="40" style="width:50;">
<label>Number of points</label>
<input id="input" type="range" min="0" max="150" value="2" style="width:240px;">
<i>Radius</i> = <output name="output" for="input">2</output>
</form>
<form>
<label for="curves">Objects</label>
<select id="curves"></select><br>
</form>
</body>
<script type="text/javascript" src="http://mbostock.github.com/d3/d3.js"></script>
<script>
/**
* Author: Brooks Mershon
* November 23, 2014
* Version 3
*
* Created as part of a project for MATH 412: Topology with Applications taught at Duke University.
*
* Copyright (c) 2014 Brooks Mershon
* The MIT License - http://opensource.org/licenses/MIT
*
*
**/
(function() {
/**
* Set up margins, SVG
* Define curves
* Define event-handlers for user inputs
*/
var margin = {top: 20, right: 20, bottom: 20, left: 20},
padding = {top: 60, right: 60, bottom: 60, left: 60},
outerWidth = 960,
outerHeight = 500,
innerWidth = outerWidth - margin.left - margin.right,
innerHeight = outerHeight - margin.top - margin.bottom,
width = innerWidth - padding.left - padding.right,
height = innerHeight - padding.top - padding.bottom;
if (self.frameElement) self.frameElement.style.height = 800 + "px";
var svg = d3.select("body").append("svg")
.attr("width", outerWidth)
.attr("height", outerHeight)
.append("g")
.attr("transform", "translate(" + margin.left + "," + margin.top + ")");
// color scale for "heatmap" indicating the time at which particular triangles are introduced using fill color
var color = d3.scale.linear()
.domain([0, 40, 60, 75, 100, 150])
.range(["#0a0", "#6c0", "#ee0", "#eb4", "#eb9", "#fff"]);
/**
* Path to sample from
**/
var epsilon_noise = 20,
n = 55;
var curves, curveNames;
d3.tsv("curves.txt", function(data) {
curves = data.map(function(d) {return d.d});
var path = svg.append("path")
.attr("d", curves[0]);
path
.attr("id", "squiggle")
.attr("stroke", "black")
.attr("stroke-opacity", 0.0)
.attr("stroke-width", 2)
.attr("fill", "none")
curveNames = data.map(function(d) {return d.name});
d3.select("#curves")
.on("change", change)
.selectAll("option")
.data(curveNames)
.enter().append("option")
.attr("value", function(d) { return d; })
.text(function(d) { return d; });
sample();
})
/*
Add legend for triangle coloring
*/
var legend = svg.selectAll(".legend")
.data(color.ticks(6).slice(1).reverse())
.enter().append("g")
.attr("class", "legend")
.attr("transform", function(d, i) { return "translate(" + (width) + "," + (100 + i * 20) + ")"; });
legend.append("rect")
.attr("width", 20)
.attr("height", 20)
.style("fill", color);
legend.append("text")
.attr("x", 26)
.attr("y", 10)
.attr("dy", ".35em")
.text(String);
svg.append("text")
.attr("class", "label")
.attr("x", width)
.attr("y", 80)
.attr("dy", ".35em")
.text("Radius");
var betti_numbers = ["?", "?"]; // data to be displayed to user with svg: text
d3.select("#input").on("input", function() { update(this.value); });
d3.select("#sample").on("click", function() { sample(); });
d3.select("#curve_checkbox").on("change", function() { checkboxUpdate(this.checked)})
d3.select("#sample-number").on("input", function() { n = this.value; sample();});
d3.select("#betti1_button").on("click", function() { updateBetti1() });
d3.select("#noise-number").on("input", function() { epsilon_noise = this.value; sample();});
var data = [],
vertices = [];
edges = [],
triangles = [],
simplices = [],
edge_values = [],
triangle_values = [];
function sample() {
d3.selectAll(".dot").remove();
d3.selectAll(".ball").remove();
var pathNode = svg.select("#squiggle").node();
var pathLength = pathNode.getTotalLength();
var gaussianNoise = d3.random.normal(0, epsilon_noise)
data = d3.range(n).map(function() {
var p = pathNode.getPointAtLength(Math.random()*pathLength)
return [p.x + gaussianNoise(), p.y + gaussianNoise()];
});
vertices = d3.range(n); // first vertex is "0"
var dot = svg.selectAll("dot")
.data(data)
.enter().append("circle")
dot
.attr("class", "dot")
.attr("cx", function(d) {return d[0]})
.attr("cy", function(d) {return d[1]})
.attr("r", 2);
var circle = svg.selectAll("ball")
.data(data)
.enter().append("circle");
circle
.attr("class", "ball")
.attr("cx", function(d) {return d[0]})
.attr("cy", function(d) {return d[1]})
.attr("r", 1);
update();
}
function indexOfPoint(array, point) {
var epsilon = 0.0000001;
var index = -1;
var length = array.length;
for(var i = 0; i < length; i++) {
var current = array[i];
if(Math.abs(point[0] - current[0]) < epsilon && Math.abs(point[1] - current[1]) < epsilon){
index = i;
break;
}
}
return index;
}
function indexOfEdge(array, edge){
var index = -1;
var length = array.length;
for(var i = 0; i < length; i++) {
var current = array[i];
if(edge[0] == current[0] && edge[1] == current[1]){
index = i;
break;
}
}
return index;
}
/*
array contains delauney triangles, where each triple is a triple of vertex indices
*/
function indexOfTriangle(array, triangle){
var index = -1;
var length = array.length;
for(var i = 0; i < length; i++) { //for each triple in delaunay triangulation
var current = array[i]; //a triple ex) [1, 4, 0], indices refer to vertices
var count = 0;
if(current.indexOf(triangle[0]) > -1 && current.indexOf(triangle[1]) > -1 && current.indexOf(triangle[2]) > -1){
index = i;
}
}
return index;
}
function update(radius) {
simplices = [];
edges = [];
triangles = [];
svg.selectAll(".edge").remove();
svg.selectAll(".triangle").remove();
var circle = d3.selectAll(".ball")
.attr("r", radius)
/*
Compute Delaunay triangulation, using this triangulation as a limiter, so that triangles are only added to the rips-complex if
they are present in the Delaunay triangulation
*/
var delaunay = d3.geom.delaunay(data); // data looks like: [[[100, 150], [200, 190], [300, 350]], [] ...]
var delaunayEdges = [];
var delaunayVertices = delaunay.slice(0);
var delaunayLength = delaunay.length;
for(var i = 0; i < delaunayLength; i++){
for(var k = 0; k < 3; k++){
var p = delaunay[i][k];
delaunay[i][k] = indexOfPoint(data, p); // replace [coordinate_x, coordinate_y, coordinate_z] with [index_x, index_y, index_z]
}
var triple = delaunay[i];
delaunayEdges.push([triple[0], triple[1]]);
delaunayEdges.push([triple[0], triple[2]]);
delaunayEdges.push([triple[1], triple[2]]);
}
/*
Check pair-wise distances, adding an edge between points if their balls encompass them AND the edge is in the Delaunay triangulation
*/
var dataLength = data.length;
for(var i = 0; i < dataLength; i++) {
for(var j = i+1; j < dataLength; j++) {
var distance = distanceBetween(data[i], data[j]);
if(radius >= distance/2 && (indexOfEdge(delaunayEdges, [i, j]) > -1 || indexOfEdge(delaunayEdges, [j, i]) > -1)) {
svg.append("line")
.attr("class", "edge")
.attr("x1", data[i][0])
.attr("x2", data[j][0])
.attr("y1", data[i][1])
.attr("y2", data[j][1]);
edges.push([[i, j], distance/2]); // add a simplex
edge_values.push(radius/2);
}
}
}
/*
Check triples of edges in array containing edges for triangles
the faces of a triangles are indices of its three edges, which will be filled in after all simplices are sorted
*/
for(var i = 0; i < dataLength; i++) {
for(var j = i+1; j < dataLength; j++) {
for(var k = j+1; k < dataLength; k++) {
var max = Math.max(distanceBetween(data[j], data[k]), Math.max(distanceBetween(data[i], data[j]), distanceBetween(data[i], data[k])));
if(indexOfTriangle(delaunay, [i, j, k]) > -1 && radius >= max/2) { //check if triple of vertices forms a Delaunay triangle
var string = data[i][0] + "," + data[i][1] + ", " + data[j][0] + "," + data[j][1] + ", " + data[k][0] + "," + data[k][1];
svg.append("polygon")
.attr("class", "triangle")
.attr("points", string)
.attr("fill", color(max/2)) // heatmap for triangles, colored to represent the radius at which they are introduced
.attr("fill-opacity", 1)
.attr("stroke", "black")
.attr("stroke-width", "0.5px")
triangles.push([[i, j, k], max/2]); // add a simplex
}
}
}
}
simplices = simplices.concat(edges).concat(triangles); // does not include vertices
/*
Sort simplices:
break ties with lower dimension simplices coming first
*/
simplices.sort(function(a, b) {
var diff = greatestEdgeLength(a) - greatestEdgeLength(b);
return (diff == 0) ? (a[0].length - b[0].length) : diff;
})
// format simplices so they are of the form:
// triangle [edge, edge, edge], edge: [vertex, vertex]
for(var i = 0; i < simplices.length; i++){
simplices[i] = simplices[i][0].slice(0);
}
//format edges
for(var i = 0; i < edges.length; i++){
edges[i] = edges[i][0].slice(0);
}
simplices = vertices.concat(simplices);
/*
replace triangle information (inclued vertices) with 1-dimensional (edge) faces
*/
var simplicesLength = simplices.length;
var verticesLength = vertices.length;
for(var i = 0; i < simplicesLength; i++) {
if(simplices[i].length == 3) { //triangle
var t = simplices[i];
// define three faces of triangles
var faces = [[t[0], t[1]], [t[0], t[2]],[t[1], t[2]]];
for(var k = 0; k < 3; k++){
for(var j = verticesLength; j < i; j++) { //start at first edge
var edge = simplices[j];
if(simplices[j].length == 2 && faces[k][0] == edge[0] && faces[k][1] == edge[1]) {
simplices[i][k] = j; // the kth face of the triangle is the simplex (an edge) at index j
}
}
}
}
}
updateBetti0();
}
function updateBetti0() {
// all vertices born at zero
var immortal_pairs = findBirthDeathPairs(d3.range(n), edges, d3.range(n).map(function() {return 0}), edge_values)[1];
betti_numbers[0] = immortal_pairs.length;
updateText();
}
function updateBetti1() {
d3.select("body").style("background-color", "#f5f5f5");
var warning = svg.append("text")
.attr("id", "warning")
.attr("x", 200)
.attr("y", 200)
.style("font-size", "50px")
.text("Computing...")
setTimeout(function(){
var info = runRipsAlgorithm();
var loops = info[0];
betti_numbers[1] = loops;
updateText();
svg.select("#warning").remove();
d3.select("body").style("background-color", "white");
}, 200);
}
function updateText() {
var text = svg.selectAll(".bettiText")
.data(betti_numbers);
text.enter().append("text")
.attr("x", function(d, i) { return width - 150; })
.attr("y", function(d, i) { return 100 + i* 60 })
.attr("dy", ".35em")
.style("font-size","34px");
var subscripts = ["₀","₁"];
text.text(function(d, i) { return "\u03B2" + subscripts[i] + " = " + d; }) // isn't that fancy?
.attr("class", "bettiText");
}
function change() {
d3.select("path").attr("d", curves[curveNames.indexOf(this.value)]);
sample();
}
function greatestEdgeLength(t) {
return t[1];
}
function distanceBetween(a, b) {
return Math.sqrt((a[0] - b[0]) * (a[0] - b[0]) + (a[1] - b[1]) * (a[1] - b[1]));
}
function checkboxUpdate(checked) {
svg.select("#squiggle")
.attr("stroke-opacity", function(d) { return (checked) ? 0.2 : 0.0});
}
/**
* finds birth/death pairs for given vertices, edges, and function values
* Union-Find
* @return {[[finite_pairs], [immortal births]]}
*/
function findBirthDeathPairs(vertices, edges, F, H) {
// edges should come in in order of their greatest endpoint f-value
edges.sort(function (a, b) {
var diff = H[a] - H[b];
});
var u = [],
pairs = [],
immortal_pairs = [];
/*
Fill vector v with all points. Each point is it's own component initially.
*/
var verticesLength = vertices.length;
for(var i = 0; i < verticesLength; i++) {
u.push(i);
}
var edgesLength = edges.length;
for(var j = 0; j < edgesLength; j++) {
var e = edges[j],
a = e[0], //initial vertex of edge
b = e[1]; //final vertex of edge
while (a != u[a]) {a = u[a]};
while (b != u[b]) {b = u[b]};
if(a < b) {
u[b] = a;
} else if (a > b) {
u[a] = b;
}
pairs.push([0, H[j]]); // [f(vertex a, f(edge)]
}
var uLength = u.length;
// Find dots of infinite persistence, tag them as dying at -1
for(var i = 0; i < uLength; i++) {
if(u[i] == i) {
immortal_pairs.push([F[i],-1, i])
}
}
return [pairs, immortal_pairs];
}
/*
Return information about 1-diagram: [num loops, [pairs]];
*/
function runRipsAlgorithm() {
var matrix = [],
edgesAndTriangles = simplices.slice(0), // includes vertices to produce square matrix - the naive implementaition
verticesLength = vertices.length;
// boundary matrix has columns of the usal matrix as each element in the array.
var boundary_matrix = [];
var edgesAndTrianglesLength = edgesAndTriangles.length;
for(var j = 0; j < edgesAndTrianglesLength; j++) { // build up columns
var column = [];
var simplex = simplices[j]; // edge or a triangle, length 2 or length 3
for(var i = 0; i < simplices.length; i++) { // all vertices, edges, triangles
if(simplex.length == 2 && (simplex[0] == i || simplex[1] == i)){ //edge contains face i (a vertex)
column.push(1);
} else if(simplex.length == 3 && (simplex[0] == i || simplex[1] == i || simplex[2] == i)) { //triangle contains face i (an edge)
column.push(1);
} else {
column.push(0);
}
}
boundary_matrix.push(column);
}
reduced = reduce(boundary_matrix);
var lowestArray = getLowestArray(reduced);
var lowestArrayLength = lowestArray.length;
var pairs = [];
var reducedLength = reduced.length;
for(var j = verticesLength; j < reducedLength; j++){
if(reduced[j].indexOf(1) < 0) { //zero-column: a cycle is born at j
var pair = [j, -1]; // default is a cycle born at simplex-j, dead at infinity
for(var k = 0; k < lowestArrayLength; k++){ //scan low-array for first lowest-'1' in j'th row
if(lowestArray[k] == j) { // cycle dead at k (first lowest 1 encountered in jth row)
pair[1] = k;
break;
}
}
pairs.push(pair);
}
}
/*
Count number of immortal pairs. This is Betti-1
*/
var numLoops = 0;
for(var i = 0; i < pairs.length; i++){
if(pairs[i][1] == -1){
numLoops++;
}
}
return [numLoops, pairs];
}
/*
matrix is a boundary 0-1 matrix where each 1 or 0 is accessed with matrix[column][row] := element
note: each elementin matrix is a column (not a row)
Return reduced data structure using lowest-1 algorithm
*/
function reduce(matrix) {
var reducedMatrix = matrix.slice(0);
var reducedMatrixLength = reducedMatrix.length;
var currentIndex = 1;
while(currentIndex < reducedMatrixLength) {
var counter = currentIndex; //if we reach 0, we have no conflicts to the left of the current column
for(var j = 0; j < currentIndex; j++){
var lowestArray = getLowestArray(reducedMatrix);
if(lowestArray[j] == lowestArray[currentIndex] && lowestArray[j] != -1){
for(var p = 0; p < reducedMatrixLength; p++){ // add two columns
reducedMatrix[currentIndex][p] = reducedMatrix[j][p] ^ reducedMatrix[currentIndex][p];
}
} else {
counter--;
}
}
if (counter == 0) {
currentIndex++;
}
}
return reducedMatrix;
}
/*
Return array containing, for each column in boundary array, the lowest-'1' in each column
*/
function getLowestArray(m){
var low = [];
for(var col = 0; col < m.length; col++) {
low.push(m[col].lastIndexOf(1));
}
return low;
}
})();
</script>
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment