Clean-up and sync with the slides
authorDmitriy Morozov <dmitriy@mrzv.org>
Sat, 16 Jun 2012 12:13:53 -0700
changeset 2 4b3728f0d920
parent 1 a0e045bf9248
child 3 ef86268a3695
Clean-up and sync with the slides
04-1-filtration.py
04-2-persistence.py
04-alpha-shapes.py
05-alpha-shapes.py
05-rips.py
06-flat-torus.py
06-rips.py
07-ls-filtration.py
08-cycle-chain.py
08-extended-persistence.py
09-alpha-shapes.py
10-circular.py
data/chain-linkages.pts
tutorial.py
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/04-1-filtration.py	Sat Jun 16 12:13:53 2012 -0700
@@ -0,0 +1,8 @@
+simplices = [([0], 1),   ([1], 2), ([0,1], 3), ([2], 4), \
+             ([1,2], 5), ([0,2], 6)]
+f = Filtration()
+for vertices, time in simplices:
+     f.append(Simplex(vertices, time))
+f.sort(dim_data_cmp)
+for s in f:
+    print s, s.data    # s.data is the time
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/04-2-persistence.py	Sat Jun 16 12:13:53 2012 -0700
@@ -0,0 +1,6 @@
+p = StaticPersistence(f)
+p.pair_simplices()
+dgms = init_diagrams(p, f)
+for i, dgm in enumerate(dgms):
+    print "Dimension:", i
+    print dgm
--- a/04-alpha-shapes.py	Tue Jun 12 22:35:53 2012 -0700
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,17 +0,0 @@
-# Alpha shapes
-points = read_points('data/trefoil.pts')
-f = Filtration()
-fill_alpha_complex(points, f)
-for s in f:
-    print s, s.data
-
-from math import sqrt
-show_complex(points, [s for s in f if sqrt(s.data[0]) < .5])
-show_complex(points, [s for s in f if sqrt(s.data[0]) < .8])
-
-f.sort(dim_data_cmp)
-p = StaticPersistence(f)
-p.pair_simplices()
-
-dgms = init_diagrams(p, f, lambda s: s.data[0])
-show_diagram(dgms)
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/05-alpha-shapes.py	Sat Jun 16 12:13:53 2012 -0700
@@ -0,0 +1,17 @@
+# Alpha shapes
+points = read_points('data/trefoil.pts')
+f = Filtration()
+fill_alpha_complex(points, f)
+for s in f:
+    print s, s.data
+
+from math import sqrt
+show_complex(points, [s for s in f if sqrt(s.data[0]) < .5])
+show_complex(points, [s for s in f if sqrt(s.data[0]) < .8])
+
+f.sort(dim_data_cmp)
+p = StaticPersistence(f)
+p.pair_simplices()
+
+dgms = init_diagrams(p, f, lambda s: sqrt(s.data[0]))
+show_diagram(dgms)
--- a/05-rips.py	Tue Jun 12 22:35:53 2012 -0700
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,19 +0,0 @@
-# Rips
-points = read_points('data/trefoil.pts')
-distances = PairwiseDistances(points)
-rips = Rips(distances)
-
-f = Filtration()
-rips.generate(2, 1.7, f.append)
-print "Number of simplices:", len(f)
-
-show_complex(points, f)
-show_complex(points, [s for s in f if rips.eval(s) < 1.6])
-
-f.sort(rips.cmp)
-
-p = StaticPersistence(f)
-p.pair_simplices()
-
-dgms = init_diagrams(p, f, rips.eval)
-show_diagram(dgms[:2])
--- a/06-flat-torus.py	Tue Jun 12 22:35:53 2012 -0700
+++ b/06-flat-torus.py	Sat Jun 16 12:13:53 2012 -0700
@@ -18,9 +18,10 @@
         sq_dist = ((x1 - x2) % 1)**2 + ((y1 - y2) % 1)**2
         return sqrt(sq_dist)
 
-rips = Rips(TorusDistances(200))
+distances = TorusDistances(250)
+rips = Rips(distances)
 f = Filtration()
-rips.generate(2, .5, f.append)
+rips.generate(2, .75, f.append)
 print "Number of simplices:", len(f)
 
 f.sort(rips.cmp)
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/06-rips.py	Sat Jun 16 12:13:53 2012 -0700
@@ -0,0 +1,19 @@
+# Rips
+points = read_points('data/trefoil.pts')
+distances = PairwiseDistances(points)
+distances = ExplicitDistances(distances)
+rips = Rips(distances)
+
+f = Filtration()
+rips.generate(2, 1.7, f.append)
+print "Number of simplices:", len(f)
+
+show_complex(points, f)
+show_complex(points, [s for s in f if rips.eval(s) < 1.6])
+
+f.sort(rips.cmp)
+p = StaticPersistence(f)
+p.pair_simplices()
+
+dgms = init_diagrams(p, f, rips.eval)
+show_diagram(dgms[:2])
--- a/07-ls-filtration.py	Tue Jun 12 22:35:53 2012 -0700
+++ b/07-ls-filtration.py	Sat Jun 16 12:13:53 2012 -0700
@@ -1,49 +1,35 @@
 # Lower-star filtration
 elephant_points, elephant_complex = read_off('data/cgal/elephant.off')
-show_complex(elephant_points, elephant_complex)
-
+# show_complex(elephant_points, elephant_complex)
 elephant_complex = closure(elephant_complex, 2)
 show_complex(elephant_points, elephant_complex)
 
 def axis(points, axis = 1):
     def value(v):
-        return points[v][1]
+        return points[v][axis]
 
     return value
 
 value = axis(elephant_points, 1)
 
 
-def max_vertex_compare(points):
+def max_vertex_compare(values):
     def max_vertex(s):
-        return max(value(v) for v in s.vertices)
+        return max(values(v) for v in s.vertices)
 
     def compare(s1, s2):
         return cmp(s1.dimension(), s2.dimension()) or cmp(max_vertex(s1), max_vertex(s2))
 
     return compare
 
-f = Filtration(elephant_complex, max_vertex_compare(elephant_points))
+f = Filtration(elephant_complex, max_vertex_compare(value))
 
-##persistence = StaticPersistence(f)
+#persistence = StaticPersistence(f)
 persistence = DynamicPersistenceChains(f)
 persistence.pair_simplices()
 
 dgms = init_diagrams(persistence, f, lambda s: max(value(v) for v in s.vertices), lambda i: i)
-
 print "Betti numbers:"
 for i,dgm in enumerate(dgms):
     print len([p for p in dgm if p[1] == float('inf')])
-
-while True:
-    pt = show_diagram(dgms)
-    if not pt: break
-    print pt
-    i = pt[2]
-    smap = persistence.make_simplex_map(f)
-    chain = [smap[ii] for ii in i.chain]
-    pair_cycle = [smap[ii] for ii in i.pair().cycle]
-    pair_chain = [smap[ii] for ii in i.pair().chain]
-    show_complex(elephant_points, subcomplex = chain)
-    show_complex(elephant_points, subcomplex = pair_cycle + pair_chain)
-    #show_complex(elephant_points, [s for s in f if s.dimension() != 1], subcomplex = cycle)
+show_diagram(dgms)
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/08-cycle-chain.py	Sat Jun 16 12:13:53 2012 -0700
@@ -0,0 +1,14 @@
+while True:
+    pt = show_diagram(dgms)
+    if not pt: break
+    print pt
+    i = pt[2]
+    smap = persistence.make_simplex_map(f)
+    chain = [smap[ii] for ii in i.chain]
+    pair_cycle = [smap[ii] for ii in i.pair().cycle]
+    pair_chain = [smap[ii] for ii in i.pair().chain]
+    show_complex(elephant_points, subcomplex = chain)
+    if pt[1] != float('inf'):
+        show_complex(elephant_points, subcomplex = pair_cycle + pair_chain)
+    #show_complex(elephant_points, [s for s in f if s.dimension() != 1], subcomplex = cycle)
+
--- a/08-extended-persistence.py	Tue Jun 12 22:35:53 2012 -0700
+++ b/08-extended-persistence.py	Sat Jun 16 12:13:53 2012 -0700
@@ -1,20 +1,27 @@
-# Extended persistence (FIXME)
+# Extended persistence
 cone = [Simplex([-1] + [v for v in s.vertices]) for s in elephant_complex]
 cone.append(Simplex([-1]))
 
-def ep_compare(points):
+def axis(points, axis = 1):
+    def value(v):
+        return points[v][axis]
+
+    return value
+value = axis(elephant_points, 1)
+
+def ep_compare(values):
     def max_vertex(s):
-        return max(points[v][2] for v in s.vertices if v != -1)
+        return max(values(v) for v in s.vertices if v != -1)
 
     def min_vertex(s):
-        return max(points[v][2] for v in s.vertices if v != -1)
+        return min(values(v) for v in s.vertices if v != -1)
 
     def compare(s1, s2):
         if s1.dimension() != s2.dimension():
             return cmp(s1.dimension(), s2.dimension())
 
         if (-1 in s1.vertices) ^ (-1 in s2.vertices):
-            return -1 if -1 in s1.vertices else 1
+            return 1 if -1 in s1.vertices else -1
         elif -1 in s1.vertices:
             return cmp(min_vertex(s1), min_vertex(s2))
         else:
@@ -23,22 +30,23 @@
     return compare
 
 f = Filtration(elephant_complex + cone)
-f.sort(ep_compare(elephant_points))
+f.sort(ep_compare(value))
 
 persistence = StaticPersistence(f)
 persistence.pair_simplices()
 
-def eval_ep(points):
+def eval_ep(values):
     def eval(s):
+        if s.dimension() == 0 and -1 in s.vertices:
+            return float('inf')
         if -1 in s.vertices:
-            return min(elephant_points[v][2] for v in s.vertices)
+            return min(values(v) for v in s.vertices if v != -1)
         else:
-            return max(elephant_points[v][2] for v in s.vertices)
-
+            return max(values(v) for v in s.vertices)
     return eval
 
-dgms = init_diagrams(persistence, f, eval_ep(elephant_points), lambda i: i)
+dgms = init_diagrams(persistence, f, eval_ep(value), lambda i: i)
 
-#show_diagram(dgms)
-print dgms
-print len(dgms)
+show_diagram(dgms)
+#print dgms
+#print len(dgms)
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/09-alpha-shapes.py	Sat Jun 16 12:13:53 2012 -0700
@@ -0,0 +1,10 @@
+from math import sqrt
+
+def alpha_shapes_diagrams(points):
+    f = Filtration()
+    fill_alpha_complex(points, f)
+    f.sort(dim_data_cmp)
+    p = StaticPersistence(f)
+    p.pair_simplices()
+    dgms = init_diagrams(p, f, lambda s: sqrt(s.data[0]))
+    return dgms
--- a/10-circular.py	Tue Jun 12 22:35:53 2012 -0700
+++ b/10-circular.py	Sat Jun 16 12:13:53 2012 -0700
@@ -9,13 +9,13 @@
 
 dgms = init_diagrams(p,f, lambda s: sqrt(s.data[0]), lambda n: n.cocycle)
 
-pt = show_diagram(dgms)
+while True:
+    pt = show_diagram(dgms)
+    if not pt: break
 
-if pt:
     #rf = Filtration((s for s in f if sqrt(s.data[0]) <= pt[0]))
     #rf = Filtration((s for s in f if sqrt(s.data[0]) < pt[1]))
     rf = Filtration((s for s in f if sqrt(s.data[0]) <= (pt[0] + pt[1])/2))
-    #rf = Filtration((s for s in f if sqrt(s.data[0]) < min(pt1[1], pt2[1], pt3[1])))
     values = circular.smooth(rf, pt[2])
     cocycle = [rf[i] for (c,i) in pt[2] if i < len(rf)]
     #print cocycle
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/data/chain-linkages.pts	Sat Jun 16 12:13:53 2012 -0700
@@ -0,0 +1,170 @@
+267.956 482.347
+257.896 493.59
+261.447 506.607
+251.387 523.176
+259.671 536.785
+263.814 536.193
+286.89 548.619
+278.015 556.312
+274.464 543.294
+299.316 547.436
+309.376 559.27
+311.743 550.986
+307.009 535.602
+299.908 517.85
+293.991 504.241
+279.79 498.323
+269.731 483.531
+268.547 475.247
+277.423 464.596
+282.157 448.619
+271.506 435.602
+275.648 418.442
+262.63 407.791
+254.938 410.158
+254.346 425.542
+242.512 428.501
+229.494 437.377
+230.677 463.412
+240.145 473.471
+238.961 456.312
+230.086 450.394
+249.612 471.104
+253.163 485.306
+256.121 394.181
+255.529 374.655
+270.322 356.903
+264.997 353.945
+251.387 359.27
+254.938 377.022
+285.115 347.436
+305.234 343.294
+314.701 348.619
+328.902 350.986
+341.92 340.335
+350.204 345.069
+358.488 346.252
+367.956 329.684
+376.24 320.809
+392.808 313.708
+399.316 298.323
+405.825 301.282
+415.885 311.341
+424.76 306.016
+441.92 304.241
+446.654 315.483
+446.654 330.276
+456.121 339.744
+460.263 358.678
+467.956 381.755
+478.606 387.673
+484.524 400.099
+470.914 391.223
+463.814 367.554
+482.157 395.957
+499.316 406.607
+515.885 413.116
+532.453 427.909
+536.003 443.886
+531.269 464.004
+530.677 481.164
+515.293 478.797
+498.725 481.755
+494.583 491.815
+486.89 496.548
+468.547 500.099
+472.689 513.708
+462.63 532.643
+459.08 541.519
+459.671 562.229
+459.08 572.88
+444.879 583.531
+429.494 599.507
+421.21 591.815
+402.275 588.264
+393.399 586.489
+380.382 572.288
+375.056 557.495
+364.997 551.578
+351.979 555.72
+344.287 546.252
+318.843 544.477
+324.76 551.578
+335.411 540.335
+457.896 520.809
+479.198 491.815
+292.808 339.152
+270.322 401.282
+292.808 394.181
+301.092 387.673
+293.399 379.389
+309.376 372.288
+309.376 364.004
+298.725 358.678
+286.299 394.773
+309.967 340.335
+305.234 326.134
+321.21 317.258
+321.21 313.708
+332.453 295.365
+347.245 304.241
+351.979 303.057
+364.997 300.099
+370.914 310.158
+372.689 314.3
+389.257 324.951
+399.316 339.152
+416.476 346.252
+438.961 342.702
+442.512 343.294
+436.595 354.536
+424.76 348.619
+408.784 345.069
+399.316 347.436
+444.287 349.211
+459.671 344.477
+474.464 340.335
+489.849 348.028
+498.133 349.803
+512.926 361.045
+515.293 375.838
+509.376 388.856
+509.967 397.732
+507.6 411.933
+497.541 421.4
+486.89 429.684
+489.257 445.661
+494.583 456.312
+497.541 468.146
+492.808 476.43
+499.908 490.631
+512.926 508.383
+515.885 521.992
+502.867 540.927
+496.358 550.394
+508.784 528.501
+509.967 506.607
+508.192 495.957
+482.748 556.312
+489.849 549.803
+464.405 553.945
+455.529 548.619
+436.595 545.661
+428.311 546.252
+412.926 555.128
+412.926 556.903
+401.092 558.087
+392.216 566.371
+392.216 566.371
+408.784 542.11
+416.476 540.927
+381.565 575.247
+370.914 590.039
+356.121 593.59
+340.145 589.448
+331.269 577.022
+321.802 570.513
+317.66 566.371
+531.861 422.584
+520.618 417.85
+534.819 447.436
--- a/tutorial.py	Tue Jun 12 22:35:53 2012 -0700
+++ b/tutorial.py	Sat Jun 16 12:13:53 2012 -0700
@@ -1,11 +1,13 @@
 execfile("01-setup.py")
 execfile("02-simplex.py")
 execfile("03-complex.py")
-execfile("04-alpha-shapes.py")
-execfile("05-rips.py")
-execfile("06-flat-torus.py")
+execfile("04-1-filtration.py")
+execfile("04-2-persistence.py")
+execfile("05-alpha-shapes.py")
+execfile("06-rips.py")
+#execfile("06-flat-torus.py")
 execfile("07-ls-filtration.py")
-execfile("08-extended-persistence.py")         # FIXME
+#execfile("08-extended-persistence.py")
 
 # TODO: Distances