1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
|
<?xml version="1.0" encoding="UTF-8"?>
<!DOCTYPE html PUBLIC "-//W3C//DTD XHTML 1.0 Strict//EN"
"http://www.w3.org/TR/xhtml1/DTD/xhtml1-strict.dtd">
<html xmlns="http://www.w3.org/1999/xhtml" xml:lang="en">
<head>
<title>GAP (HAP commands) - Chapter 5: Topological data analysis</title>
<meta http-equiv="content-type" content="text/html; charset=UTF-8" />
<meta name="generator" content="GAPDoc2HTML" />
<link rel="stylesheet" type="text/css" href="manual.css" />
<script src="manual.js" type="text/javascript"></script>
<script type="text/javascript">overwriteStyle();</script>
</head>
<body class="chap5" onload="jscontent()">
<div class="chlinktop"><span class="chlink1">Goto Chapter: </span><a href="chap0.html">Top</a> <a href="chap1.html">1</a> <a href="chap2.html">2</a> <a href="chap3.html">3</a> <a href="chap4.html">4</a> <a href="chap5.html">5</a> <a href="chap6.html">6</a> <a href="chap7.html">7</a> <a href="chap8.html">8</a> <a href="chap9.html">9</a> <a href="chap10.html">10</a> <a href="chap11.html">11</a> <a href="chap12.html">12</a> <a href="chap13.html">13</a> <a href="chap14.html">14</a> <a href="chap15.html">15</a> <a href="chap16.html">16</a> <a href="chapBib.html">Bib</a> <a href="chapInd.html">Ind</a> </div>
<div class="chlinkprevnexttop"> <a href="chap0.html">[Top of Book]</a> <a href="chap0.html#contents">[Contents]</a> <a href="chap4.html">[Previous Chapter]</a> <a href="chap6.html">[Next Chapter]</a> </div>
<p id="mathjaxlink" class="pcenter"><a href="chap5_mj.html">[MathJax on]</a></p>
<p><a id="X7B7E077887694A9F" name="X7B7E077887694A9F"></a></p>
<div class="ChapSects"><a href="chap5.html#X7B7E077887694A9F">5 <span class="Heading">Topological data analysis</span></a>
<div class="ContSect"><span class="tocline"><span class="nocss"> </span><a href="chap5.html#X80A70B20873378E0">5.1 <span class="Heading">Persistent homology </span></a>
</span>
<div class="ContSSBlock">
<span class="ContSS"><br /><span class="nocss"> </span><a href="chap5.html#X7D512DA37F789B4C">5.1-1 <span class="Heading">Background to the data</span></a>
</span>
</div></div>
<div class="ContSect"><span class="tocline"><span class="nocss"> </span><a href="chap5.html#X849556107A23FF7B">5.2 <span class="Heading">Mapper clustering</span></a>
</span>
<div class="ContSSBlock">
<span class="ContSS"><br /><span class="nocss"> </span><a href="chap5.html#X7D512DA37F789B4C">5.2-1 <span class="Heading">Background to the data</span></a>
</span>
</div></div>
<div class="ContSect"><span class="tocline"><span class="nocss"> </span><a href="chap5.html#X7BBDE0567DB8C5DA">5.3 <span class="Heading">Some tools for handling pure complexes</span></a>
</span>
</div>
<div class="ContSect"><span class="tocline"><span class="nocss"> </span><a href="chap5.html#X79616D12822FDB9A">5.4 <span class="Heading">Digital image analysis and persistent homology</span></a>
</span>
<div class="ContSSBlock">
<span class="ContSS"><br /><span class="nocss"> </span><a href="chap5.html#X8066F9B17B78418E">5.4-1 <span class="Heading">Naive example of image segmentation by automatic thresholding</span></a>
</span>
<span class="ContSS"><br /><span class="nocss"> </span><a href="chap5.html#X7E6436E0856761F2">5.4-2 <span class="Heading">Refining the filtration</span></a>
</span>
<span class="ContSS"><br /><span class="nocss"> </span><a href="chap5.html#X7D512DA37F789B4C">5.4-3 <span class="Heading">Background to the data</span></a>
</span>
</div></div>
<div class="ContSect"><span class="tocline"><span class="nocss"> </span><a href="chap5.html#X7A8224DA7B00E0D9">5.5 <span class="Heading">A second example of digital image segmentation</span></a>
</span>
</div>
<div class="ContSect"><span class="tocline"><span class="nocss"> </span><a href="chap5.html#X8290E7D287F69B98">5.6 <span class="Heading">A third example of digital image segmentation</span></a>
</span>
</div>
<div class="ContSect"><span class="tocline"><span class="nocss"> </span><a href="chap5.html#X7957F329835373E9">5.7 <span class="Heading">Naive example of digital image contour extraction</span></a>
</span>
</div>
<div class="ContSect"><span class="tocline"><span class="nocss"> </span><a href="chap5.html#X7D2CC9CB85DF1BAF">5.8 <span class="Heading">Alternative approaches to computing persistent homology</span></a>
</span>
<div class="ContSSBlock">
<span class="ContSS"><br /><span class="nocss"> </span><a href="chap5.html#X86FD0A867EC9E64F">5.8-1 <span class="Heading">Non-trivial cup product</span></a>
</span>
<span class="ContSS"><br /><span class="nocss"> </span><a href="chap5.html#X783EF0F17B629C46">5.8-2 <span class="Heading">Explicit homology generators</span></a>
</span>
</div></div>
<div class="ContSect"><span class="tocline"><span class="nocss"> </span><a href="chap5.html#X80D0D8EB7BCD05E9">5.9 <span class="Heading">Knotted proteins</span></a>
</span>
</div>
<div class="ContSect"><span class="tocline"><span class="nocss"> </span><a href="chap5.html#X87AF06677F05C624">5.10 <span class="Heading">Random simplicial complexes</span></a>
</span>
</div>
<div class="ContSect"><span class="tocline"><span class="nocss"> </span><a href="chap5.html#X875EE92F7DBA1E27">5.11 <span class="Heading">Computing homology of a clique complex (Vietoris-Rips complex) </span></a>
</span>
</div>
</div>
<h3>5 <span class="Heading">Topological data analysis</span></h3>
<p><a id="X80A70B20873378E0" name="X80A70B20873378E0"></a></p>
<h4>5.1 <span class="Heading">Persistent homology </span></h4>
<p>Pairwise distances between <span class="SimpleMath">74</span> points from some metric space have been recorded and stored in a <span class="SimpleMath">74× 74</span> matrix <span class="SimpleMath">D</span>. The following commands load the matrix, construct a filtration of length <span class="SimpleMath">100</span> on the first two dimensions of the assotiated clique complex (also known as the <em>Vietoris-Rips Complex</em>), and display the resulting degree <span class="SimpleMath">0</span> persistent homology as a barcode. A single bar with label <span class="SimpleMath">n</span> denotes <span class="SimpleMath">n</span> bars with common starting point and common end point.</p>
<div class="example"><pre>
<span class="GAPprompt">gap></span> <span class="GAPinput">file:=HapFile("data253a.txt");;</span>
<span class="GAPprompt">gap></span> <span class="GAPinput">Read(file);</span>
<span class="GAPprompt">gap></span> <span class="GAPinput">G:=SymmetricMatrixToFilteredGraph(D,100);</span>
Filtered graph on 74 vertices.
<span class="GAPprompt">gap></span> <span class="GAPinput">K:=FilteredRegularCWComplex(CliqueComplex(G,2));</span>
Filtered regular CW-complex of dimension 2
<span class="GAPprompt">gap></span> <span class="GAPinput">P:=PersistentBettiNumbers(K,0);;</span>
<span class="GAPprompt">gap></span> <span class="GAPinput">BarCodeCompactDisplay(P);</span>
</pre></div>
<p><img src="images/bar0.png" align="center" height="60" alt="degree 0 barcode"/></p>
<p>The first 54 terms in the filtration each have 74 path components -- one for each point in the sample. During the next 9 filtration terms the number of path components reduces, meaning that sample points begin to coalesce due to the formation of edges in the simplicial complexes. Then, two path components persist over an interval of 18 filtration terms, before they eventually coalesce.</p>
<p>The next commands display the resulting degree <span class="SimpleMath">1</span> persistent homology as a barcode.</p>
<div class="example"><pre>
<span class="GAPprompt">gap></span> <span class="GAPinput">P:=PersistentBettiNumbers(K,1);;</span>
<span class="GAPprompt">gap></span> <span class="GAPinput">BarCodeCompactDisplay(P);</span>
</pre></div>
<p><img src="images/bar1.png" align="center" height="120" alt="degree 1 bar code"/></p>
<p>Interpreting short bars as noise, we see for instance that the <span class="SimpleMath">65</span>th term in the filtration could be regarded as noiseless and belonging to a "stable interval" in the filtration with regards to first and second homology functors. The following command displays (up to homotopy) the <span class="SimpleMath">1</span> skeleton of the simplicial complex arizing as the <span class="SimpleMath">65</span>-th term in the filtration on the clique complex.</p>
<div class="example"><pre>
<span class="GAPprompt">gap></span> <span class="GAPinput">Y:=FiltrationTerm(K,65);</span>
Regular CW-complex of dimension 1
<span class="GAPprompt">gap></span> <span class="GAPinput">Display(HomotopyGraph(Y));</span>
</pre></div>
<p><img src="images/twocircles.png" align="center" height="300" alt="1-skeleton"/></p>
<p>These computations suggest that the dataset contains two persistent path components (or clusters), and that each path component is in some sense periodic. The final command displays one possible representation of the data as points on two circles.</p>
<p><a id="X7D512DA37F789B4C" name="X7D512DA37F789B4C"></a></p>
<h5>5.1-1 <span class="Heading">Background to the data</span></h5>
<p>Each point in the dataset was an image consisting of <span class="SimpleMath">732× 761</span> pixels. This point was regarded as a vector in <span class="SimpleMath">R^557052= R^732× 761</span> and the matrix <span class="SimpleMath">D</span> was constructed using the Euclidean metric. The images were the following:</p>
<p><img src="images/letters.png" align="center" height="220" alt="letters"/></p>
<p><a id="X849556107A23FF7B" name="X849556107A23FF7B"></a></p>
<h4>5.2 <span class="Heading">Mapper clustering</span></h4>
<p>The following example reads in a set <span class="SimpleMath">S</span> of vectors of rational numbers. It uses the Euclidean distance <span class="SimpleMath">d(u,v)</span> between vectors. It fixes some vector <span class="SimpleMath">u_0∈ S</span> and uses the associated function <span class="SimpleMath">f: D→ [0,b] ⊂ R, v↦ d(u_0,v)</span>. In addition, it uses an open cover of the interval <span class="SimpleMath">[0,b]</span> consisting of <span class="SimpleMath">100</span> uniformly distributed overlapping open subintervals of radius <span class="SimpleMath">r=29</span>. It also uses a simple clustering algorithm implemented in the function <code class="code">cluster</code>.</p>
<p>These ingredients are input into the Mapper clustering procedure to produce a simplicial complex <span class="SimpleMath">M</span> which is intended to be a representation of the data. The complex <span class="SimpleMath">M</span> is <span class="SimpleMath">1</span>-dimensional and the final command uses GraphViz software to visualize the graph. The nodes of this simplicial complex are "buckets" containing data points. A data point may reside in several buckets. The number of points in the bucket determines the size of the node. Two nodes are connected by an edge when they contain common data points.</p>
<div class="example"><pre>
<span class="GAPprompt">gap></span> <span class="GAPinput">file:=HapFile("data134.txt");;</span>
<span class="GAPprompt">gap></span> <span class="GAPinput">Read(file);</span>
<span class="GAPprompt">gap></span> <span class="GAPinput">dx:=EuclideanApproximatedMetric;;</span>
<span class="GAPprompt">gap></span> <span class="GAPinput">dz:=EuclideanApproximatedMetric;;</span>
<span class="GAPprompt">gap></span> <span class="GAPinput">L:=List(S,x->Maximum(List(S,y->dx(x,y))));;</span>
<span class="GAPprompt">gap></span> <span class="GAPinput">n:=Position(L,Minimum(L));;</span>
<span class="GAPprompt">gap></span> <span class="GAPinput">f:=function(x); return [dx(S[n],x)]; end;;</span>
<span class="GAPprompt">gap></span> <span class="GAPinput">P:=30*[0..100];; P:=List(P, i->[i]);;</span>
<span class="GAPprompt">gap></span> <span class="GAPinput">r:=29;;</span>
<span class="GAPprompt">gap></span> <span class="GAPinput">epsilon:=75;;</span>
<span class="GAPprompt">gap></span> <span class="GAPinput"> cluster:=function(S)</span>
<span class="GAPprompt">></span> <span class="GAPinput"> local Y, P, C;</span>
<span class="GAPprompt">></span> <span class="GAPinput"> if Length(S)=0 then return S; fi;</span>
<span class="GAPprompt">></span> <span class="GAPinput"> Y:=VectorsToOneSkeleton(S,epsilon,dx);</span>
<span class="GAPprompt">></span> <span class="GAPinput"> P:=PiZero(Y);</span>
<span class="GAPprompt">></span> <span class="GAPinput"> C:=Classify([1..Length(S)],P[2]);</span>
<span class="GAPprompt">></span> <span class="GAPinput"> return List(C,x->S{x});</span>
<span class="GAPprompt">></span> <span class="GAPinput"> end;;</span>
<span class="GAPprompt">gap></span> <span class="GAPinput">M:=Mapper(S,dx,f,dz,P,r,cluster);</span>
Simplicial complex of dimension 1.
<span class="GAPprompt">gap></span> <span class="GAPinput">Display(GraphOfSimplicialComplex(M));</span>
</pre></div>
<p><img src="images/mapper.png" align="center" height="300" alt="Mapper graph"/></p>
<p><a id="X7D512DA37F789B4C" name="X7D512DA37F789B4C"></a></p>
<h5>5.2-1 <span class="Heading">Background to the data</span></h5>
<p>The datacloud <span class="SimpleMath">S</span> consists of the <span class="SimpleMath">400</span> points in the plane shown in the following picture.</p>
<p><img src="images/mappercloud.png" align="center" height="400" alt="data cloud"/></p>
<p><a id="X7BBDE0567DB8C5DA" name="X7BBDE0567DB8C5DA"></a></p>
<h4>5.3 <span class="Heading">Some tools for handling pure complexes</span></h4>
<p>A CW-complex <span class="SimpleMath">X</span> is said to be <em>pure</em> if all of its top-dimensional cells have a common dimension. There are instances where such a space <span class="SimpleMath">X</span> provides a convenient ambient space whose subspaces can be used to model experimental data. For instance, the plane <span class="SimpleMath">X= R^2</span> admits a pure regular CW-structure whose <span class="SimpleMath">2</span>-cells are open unit squares with integer coordinate vertices. An alternative, and sometimes preferrable, pure regular CW-structure on <span class="SimpleMath">R^2</span> is one where the <span class="SimpleMath">2</span>-cells are all reguar hexagons with sides of unit length. Any digital image can be thresholded to produce a black-white image and this black-white image can naturally be regared as a finite pure cellular subcomplex of either of the two proposed CW-structures on <span class="SimpleMath">R^2</span>. Analogously, thresholding can be used to represent <span class="SimpleMath">3</span>-dimensional greyscale images as finite pure cellular subspaces of cubical or permutahedral CW-structures on <span class="SimpleMath">R^3</span>, and to represent RGB colour photographs as analogous subcomplexes of <span class="SimpleMath">R^5</span>.</p>
<p>In this section we list a few functions for performing basic operations on <span class="SimpleMath">n</span>-dimensional pure cubical and pure permutahedral finite subcomplexes <span class="SimpleMath">M</span> of <span class="SimpleMath">X=R^n</span>. We refer to <span class="SimpleMath">M</span> simply as a <em>pure complex</em>. In subsequent sections we demonstrate how these few functions on pure complexes allow for in-depth analysis of experimental data.</p>
<p>(<strong class="button">Aside.</strong> The basic operations could equally well be implemented for other CW-decompositions of <span class="SimpleMath">X= R^n</span> such as the regular CW-decompositions arising as the tessellations by a fundamental domain of a Bieberbach group (=torsion free crytallographic group). Moreover, the basic operations could also be implemented for other manifolds such as an <span class="SimpleMath">n</span>-torus <span class="SimpleMath">X=S^1× S^1 × ⋯ × S^1</span> or <span class="SimpleMath">n</span>-sphere <span class="SimpleMath">X=S^n</span> or for <span class="SimpleMath">X</span> the universal cover of some interesting hyperbolic <span class="SimpleMath">3</span>-manifold. An example use of the ambient manifold <span class="SimpleMath">X=S^1× S^1× S^1</span> could be for the construction of a cellular subspace recording the time of day, day of week and week of the year of crimes committed in a population.)</p>
<p><strong class="button">Basic operations returning pure complexes.</strong> ( Function descriptions available <span class="URL"><a href="../doc/chap1_mj.html#X7FD50DF6782F00A0">here</a></span>.)</p>
<ul>
<li><p><code class="code">PureCubicalComplex(binary array)</code></p>
</li>
<li><p><code class="code">PurePermutahedralComplex(binary array)</code></p>
</li>
<li><p><code class="code">ReadImageAsPureCubicalComplex(file,threshold)</code></p>
</li>
<li><p><code class="code">ReadImageSquenceAsPureCubicalComplex(file,threshold)</code></p>
</li>
<li><p><code class="code">PureComplexBoundary(M)</code></p>
</li>
<li><p><code class="code">PureComplexComplement(M)</code></p>
</li>
<li><p><code class="code">PureComplexRandomCell(M)</code></p>
</li>
<li><p><code class="code">PureComplexThickened(M)</code></p>
</li>
<li><p><code class="code">ContractedComplex(M, optional subcomplex of M)</code></p>
</li>
<li><p><code class="code">ExpandedComplex(M, optional supercomplex of M)</code></p>
</li>
<li><p><code class="code">PureComplexUnion(M,N)</code></p>
</li>
<li><p><code class="code">PureComplexIntersection(M,N)</code></p>
</li>
<li><p><code class="code">PureComplexDifference(M,N)</code></p>
</li>
<li><p><code class="code">FiltrationTerm(F,n)</code></p>
</li>
</ul>
<p><strong class="button">Basic operations returning filtered pure complexes.</strong></p>
<ul>
<li><p><code class="code">PureComplexThickeningFiltration(M,length)</code></p>
</li>
<li><p><code class="code">ReadImageAsFilteredPureCubicalComplex(file,length)</code></p>
</li>
</ul>
<p><a id="X79616D12822FDB9A" name="X79616D12822FDB9A"></a></p>
<h4>5.4 <span class="Heading">Digital image analysis and persistent homology</span></h4>
<p>The following example reads in a digital image as a filtered pure cubical complexex. The filtration is obtained by thresholding at a sequence of uniformly spaced values on the greyscale range. The persistent homology of this filtered complex is calculated in degrees <span class="SimpleMath">0</span> and <span class="SimpleMath">1</span> and displayed as two barcodes.</p>
<div class="example"><pre>
<span class="GAPprompt">gap></span> <span class="GAPinput">file:=HapFile("image1.3.2.png");;</span>
<span class="GAPprompt">gap></span> <span class="GAPinput">F:=ReadImageAsFilteredPureCubicalComplex(file,40);</span>
Filtered pure cubical complex of dimension 2.
<span class="GAPprompt">gap></span> <span class="GAPinput">P:=PersistentBettiNumbers(F,0);;</span>
<span class="GAPprompt">gap></span> <span class="GAPinput">BarCodeCompactDisplay(P);</span>
</pre></div>
<p><img src="images/imbar0.gif" align="center" height="400" alt="barcode"/></p>
<div class="example"><pre>
<span class="GAPprompt">gap></span> <span class="GAPinput">P:=PersistentBettiNumbers(F,1);;</span>
<span class="GAPprompt">gap></span> <span class="GAPinput">BarCodeCompactDisplay(P);</span>
</pre></div>
<p><img src="images/imbar1.gif" align="center" height="400" alt="barcode"/></p>
<p>The <span class="SimpleMath">20</span> persistent bars in the degree <span class="SimpleMath">0</span> barcode suggest that the image has <span class="SimpleMath">20</span> objects. The degree <span class="SimpleMath">1</span> barcode suggests that there are <span class="SimpleMath">14</span> (or possibly <span class="SimpleMath">17</span>) holes in these <span class="SimpleMath">20</span> objects.</p>
<p><a id="X8066F9B17B78418E" name="X8066F9B17B78418E"></a></p>
<h5>5.4-1 <span class="Heading">Naive example of image segmentation by automatic thresholding</span></h5>
<p>Assuming that short bars and isolated points in the barcodes represent noise while long bars represent essential features, a "noiseless" representation of the image should correspond to a term in the filtration corresponding to a column in the barcode incident with all the long bars but incident with no short bars or isolated points. There is no noiseless term in the above filtration of length 40. However (in conjunction with the next subsection) the following commands confirm that the 64th term in the filtration of length 500 is such a term and display this term as a binary image.</p>
<div class="example"><pre>
<span class="GAPprompt">gap></span> <span class="GAPinput">F:=ReadImageAsFilteredPureCubicalComplex(file,500);;</span>
<span class="GAPprompt">gap></span> <span class="GAPinput">Y:=FiltrationTerm(F,64); </span>
Pure cubical complex of dimension 2.
<span class="GAPprompt">gap></span> <span class="GAPinput">BettiNumber(Y,0);</span>
20
<span class="GAPprompt">gap></span> <span class="GAPinput">BettiNumber(Y,1);</span>
14
<span class="GAPprompt">gap></span> <span class="GAPinput">Display(Y);</span>
</pre></div>
<p><img src="images/binaryimage.png" align="center" height="400" alt="binary image"/></p>
<p><a id="X7E6436E0856761F2" name="X7E6436E0856761F2"></a></p>
<h5>5.4-2 <span class="Heading">Refining the filtration</span></h5>
<p>The first filtration for the image has 40 terms. One may wish to investigate a filtration with more terms, say 500 terms, with a view to analysing, say, those 1-cycles that are born by term 25 of the filtration and that die between terms 50 and 60. The following commands produce the relevant barcode showing that there is precisely one such 1-cycle.</p>
<div class="example"><pre>
<span class="GAPprompt">gap></span> <span class="GAPinput">F:=ReadImageAsFilteredPureCubicalComplex(file,500);;</span>
<span class="GAPprompt">gap></span> <span class="GAPinput">L:=[20,60,61,62,63,64,65,66,67,68,69,70];; </span>
<span class="GAPprompt">gap></span> <span class="GAPinput">T:=FiltrationTerms(F,L);;</span>
<span class="GAPprompt">gap></span> <span class="GAPinput">P0:=PersistentBettiNumbers(T,0);;</span>
<span class="GAPprompt">gap></span> <span class="GAPinput">BarCodeCompactDisplay(P0);</span>
<span class="GAPprompt">gap></span> <span class="GAPinput">P1:=PersistentBettiNumbers(T,1);;</span>
<span class="GAPprompt">gap></span> <span class="GAPinput">BarCodeCompactDisplay(P1);</span>
</pre></div>
<p><span class="SimpleMath">β_0</span>:</p>
<p><img src="images/refinedbc0.gif" align="center" height="100" alt="bar code"/></p>
<p><span class="SimpleMath">β_1</span>:</p>
<p><img src="images/refinedbc.gif" align="center" height="200" alt="bar code"/></p>
<p><a id="X7D512DA37F789B4C" name="X7D512DA37F789B4C"></a></p>
<h5>5.4-3 <span class="Heading">Background to the data</span></h5>
<p>The following image was used in the example.</p>
<p><img src="../tst/examples/image1.3.2.png" align="center" height="400" alt="barcode"/></p>
<p><a id="X7A8224DA7B00E0D9" name="X7A8224DA7B00E0D9"></a></p>
<h4>5.5 <span class="Heading">A second example of digital image segmentation</span></h4>
<p>In order to automatically count the number of coins in this picture</p>
<p><img src="images/my_coins.png" align="center" height="400" alt="collection of coins"/></p>
<p>we can load the image as a filtered pure cubical complex <span class="SimpleMath">F</span> of filtration length 40 say, and observe the degree zero persistent Betti numbers to establish that the 28-th term or so of <span class="SimpleMath">F</span> seems to be 'noise free' in degree zero. We can then set <span class="SimpleMath">M</span> equal to the 28-th term of <span class="SimpleMath">F</span> and thicken <span class="SimpleMath">M</span> a couple of times say to remove any tiny holes it may have. We can then construct the complement <span class="SimpleMath">C</span> of <span class="SimpleMath">M</span>. Then we can construct a 'neighbourhood thickening' filtration <span class="SimpleMath">T</span> of <span class="SimpleMath">C</span> with say <span class="SimpleMath">50</span> consecutive thickenings. The degree one persistent barcode for <span class="SimpleMath">T</span> has <span class="SimpleMath">24</span> long bars, suggesting that the original picture consists of <span class="SimpleMath">24</span> coins.</p>
<div class="example"><pre>
<span class="GAPprompt">gap></span> <span class="GAPinput">F:=ReadImageAsFilteredPureCubicalComplex("my_coins.png",40);;</span>
<span class="GAPprompt">gap></span> <span class="GAPinput">M:=FiltrationTerm(F,24);; #Chosen after viewing degree 0 barcode for F</span>
<span class="GAPprompt">gap></span> <span class="GAPinput">M:=PureComplexThickened(M);;</span>
<span class="GAPprompt">gap></span> <span class="GAPinput">M:=PureComplexThickened(M);;</span>
<span class="GAPprompt">gap></span> <span class="GAPinput">C:=PureComplexComplement(M);;</span>
<span class="GAPprompt">gap></span> <span class="GAPinput">T:=ThickeningFiltration(C,50);;</span>
<span class="GAPprompt">gap></span> <span class="GAPinput">P:=PersistentBettiNumbers(T,1);;</span>
<span class="GAPprompt">gap></span> <span class="GAPinput">BarCodeCompactDisplay(P);</span>
</pre></div>
<p><img src="images/coinsbettizero.gif" align="center" height="400" alt="barcode"/></p>
<p>The pure cubical complex <code class="code">W:=PureComplexComplement(FiltrationTerm(T,25))</code> has the correct number of path components, namely <span class="SimpleMath">25</span>, but its path components are very much subsets of the regions in the image corresponding to coins. The complex <span class="SimpleMath">W</span> can be thickened repeatedly, subject to no two path components being allowed to merge, in order to obtain a more realistic image segmentation with path components corresponding more closely to coins. This is done in the follow commands which use a makeshift function <code class="code">Basins(L)</code> available <span class="URL"><a href="tutex/basins.g">here</a></span>. The commands essentially implement a standard watershed segmentation algorithm but do so by using the language of filtered pure cubical complexes.</p>
<div class="example"><pre>
<span class="GAPprompt">gap></span> <span class="GAPinput">W:=PureComplexComplement(FiltrationTerm(T,25));;</span>
<span class="GAPprompt">gap></span> <span class="GAPinput">L:=[];;</span>
<span class="GAPprompt">gap></span> <span class="GAPinput">for i in [1..PathComponentOfPureComplex(W,0)] do</span>
<span class="GAPprompt">gap></span> <span class="GAPinput">P:=PathComponentOfPureComplex(W,i);;</span>
<span class="GAPprompt">gap></span> <span class="GAPinput">Q:=ThickeningFiltration(P,50,M);;</span>
<span class="GAPprompt">gap></span> <span class="GAPinput">Add(L,Q);;</span>
<span class="GAPprompt">gap></span> <span class="GAPinput">od;;</span>
<span class="GAPprompt">gap></span> <span class="GAPinput">B:=Basins(L);</span>
<span class="GAPprompt">gap></span> <span class="GAPinput">Display(B);</span>
</pre></div>
<p><img src="images/segmented_coins.png" align="center" height="400" alt="segmented coins"/></p>
<p><a id="X8290E7D287F69B98" name="X8290E7D287F69B98"></a></p>
<h4>5.6 <span class="Heading">A third example of digital image segmentation</span></h4>
<p>The following image is number 3096 in the <span class="URL"><a href="https://www2.eecs.berkeley.edu/Research/Projects/CS/vision/bsds/">BSDS500 database of images</a></span> <a href="chapBib.html#biBMartinFTM01">[MFTM01]</a>.</p>
<p><img src="images/3096.jpg" align="center" height="200" alt="image 3096 from BSDS500"/></p>
<p>A common first step in segmenting such an image is to appropriately threshold the corresponding gradient image.</p>
<p><img src="images/3096b.jpg" align="center" height="200" alt="gradient image"/> <img src="images/3096points.png" align="center" height="200" alt="thresholded gradient image"/></p>
<p>The following commands use the thresholded gradient image to produce an outline of the aeroplane. The outline is a pure cubical complex with one path component and with first Betti number equal to 1.</p>
<div class="example"><pre>
<span class="GAPprompt">gap></span> <span class="GAPinput">file:=Filename(DirectoriesPackageLibrary("HAP"),"../tutorial/images/3096b.jpg");;</span>
<span class="GAPprompt">gap></span> <span class="GAPinput">F:=ReadImageAsFilteredPureCubicalComplex(file,30);;</span>
<span class="GAPprompt">gap></span> <span class="GAPinput">F:=ComplementOfFilteredPureCubicalComplex(F);;</span>
<span class="GAPprompt">gap></span> <span class="GAPinput">M:=FiltrationTerm(F,27);; #Thickening chosen based on degree 0 barcode</span>
<span class="GAPprompt">gap></span> <span class="GAPinput">Display(M);;</span>
<span class="GAPprompt">gap></span> <span class="GAPinput">P:=List([1..BettiNumber(M,0)],n->PathComponentOfPureComplex(M,n));;</span>
<span class="GAPprompt">gap></span> <span class="GAPinput">P:=Filtered(P,m->Size(m)>10);;</span>
<span class="GAPprompt">gap></span> <span class="GAPinput">M:=P[1];;</span>
<span class="GAPprompt">gap></span> <span class="GAPinput">for m in P do</span>
<span class="GAPprompt">></span> <span class="GAPinput">M:=PureComplexUnion(M,m);;</span>
<span class="GAPprompt">></span> <span class="GAPinput">od;</span>
<span class="GAPprompt">gap></span> <span class="GAPinput">T:=ThickeningFiltration(M,50);;</span>
<span class="GAPprompt">gap></span> <span class="GAPinput">BettiNumber(FiltrationTerm(T,11),0);</span>
1
<span class="GAPprompt">gap></span> <span class="GAPinput">BettiNumber(FiltrationTerm(T,11),1);</span>
1
<span class="GAPprompt">gap></span> <span class="GAPinput">BettiNumber(FiltrationTerm(T,12),1);</span>
0
<span class="GAPprompt">gap></span> <span class="GAPinput">#Confirmation that 11-th filtration term has one hole and the 12-th term is contractible.</span>
<span class="GAPprompt">gap></span> <span class="GAPinput">C:=FiltrationTerm(T,11);;</span>
<span class="GAPprompt">gap></span> <span class="GAPinput">for n in Reversed([1..10]) do</span>
<span class="GAPprompt">></span> <span class="GAPinput">C:=ContractedComplex(C,FiltrationTerm(T,n));</span>
<span class="GAPprompt">></span> <span class="GAPinput">od;</span>
<span class="GAPprompt">gap></span> <span class="GAPinput">C:=PureComplexBoundary(PureComplexThickened(C));;</span>
<span class="GAPprompt">gap></span> <span class="GAPinput">H:=HomotopyEquivalentMinimalPureCubicalSubcomplex(FiltrationTerm(T,12),C);;</span>
<span class="GAPprompt">gap></span> <span class="GAPinput">B:=ContractedComplex(PureComplexBoundary(H));;</span>
<span class="GAPprompt">gap></span> <span class="GAPinput">Display(B);</span>
</pre></div>
<p><img src="images/3096final.png" align="center" height="200" alt="outline of aeroplane"/></p>
<p><a id="X7957F329835373E9" name="X7957F329835373E9"></a></p>
<h4>5.7 <span class="Heading">Naive example of digital image contour extraction</span></h4>
<p>The following greyscale image is available from the <span class="URL"><a href="http://www.ipol.im/pub/art/2014/74/FrechetAndConnectedCompDemo.tgz">online appendix</a></span> to the paper <a href="chapBib.html#biBcoeurjolly">[CKL14]</a>.</p>
<p><img src="images/circularGradient.png" align="center" height="250" alt="circular gradient image"/></p>
<p>The following commands produce a picture of contours from this image based on greyscale values. They also produce a picture of just the closed contours (the non-closed contours having been homotopy collapsed to points).</p>
<div class="example"><pre>
<span class="GAPprompt">gap></span> <span class="GAPinput">file:=Filename(DirectoriesPackageLibrary("HAP"),"../tutorial/images/circularGradient.png");;</span>
<span class="GAPprompt">gap></span> <span class="GAPinput">L:=[];; </span>
<span class="GAPprompt">gap></span> <span class="GAPinput">for n in [1..15] do</span>
<span class="GAPprompt">></span> <span class="GAPinput">M:=ReadImageAsPureCubicalComplex(file,n*30000);</span>
<span class="GAPprompt">></span> <span class="GAPinput">M:=PureComplexBoundary(M);;</span>
<span class="GAPprompt">></span> <span class="GAPinput">Add(L,M);</span>
<span class="GAPprompt">></span> <span class="GAPinput">od;;</span>
<span class="GAPprompt">gap></span> <span class="GAPinput">C:=L[1];;</span>
<span class="GAPprompt">gap></span> <span class="GAPinput">for n in [2..Length(L)] do C:=PureComplexUnion(C,L[n]); od;</span>
<span class="GAPprompt">gap></span> <span class="GAPinput">Display(C);</span>
<span class="GAPprompt">gap></span> <span class="GAPinput">Display(ContractedComplex(C));</span>
</pre></div>
<p>Contours from the above greyscale image:</p>
<p><img src="images/contours.png" align="center" height="250" alt="contours image"/></p>
<p>Closed contours from the above greyscale image:</p>
<p><img src="images/closedcontours.png" align="center" height="250" alt="closedcontours image"/></p>
<p>Very similar results are obtained when applied to the file <code class="code">circularGradientNoise.png</code>, containing noise, available from the <span class="URL"><a href="http://www.ipol.im/pub/art/2014/74/FrechetAndConnectedCompDemo.tgz">online appendix</a></span> to the paper <a href="chapBib.html#biBcoeurjolly">[CKL14]</a>.</p>
<p>The number of distinct "light sources" in the image can be read from the countours. Alternatively, this number can be read directly from the barcode produced by the following commands.</p>
<div class="example"><pre>
<span class="GAPprompt">gap></span> <span class="GAPinput">F:=ReadImageAsFilteredPureCubicalComplex(file,20);;</span>
<span class="GAPprompt">gap></span> <span class="GAPinput">P:=PersistentBettiNumbersAlt(F,1);;</span>
<span class="GAPprompt">gap></span> <span class="GAPinput">BarCodeCompactDisplay(P);</span>
</pre></div>
<p><img src="images/bccircularGradient.png" align="center" height="250" alt="closedcontours image"/></p>
<p>The seventeen bars in the barcode correspond to seventeen light sources. The length of a bar is a measure of the "persistence" of the corresponding light source. A long bar may initially represent a cluster of several lights whose members may eventually be distinguished from each other as new bars (or persistent homology classes) are created.</p>
<p>Here the command <code class="code">PersistentBettiNumbersAlt</code> has been used. This command is explained in the following section.</p>
<p>The follwowing commands use a watershed method to partition the digital image into regions, one region per light source. A makeshift function <code class="code">Basins(L)</code>, available <span class="URL"><a href="tutex/basins.g">here</a></span>, is called. (The efficiency of the example could be easily improved. For simplicity it uses generic commands which, in principle, can be applied to cubical or permutarhedral complexes of higher dimensions.)</p>
<div class="example"><pre>
<span class="GAPprompt">gap></span> <span class="GAPinput">file:=Filename(DirectoriesPackageLibrary("HAP"),"../tutorial/images/circularGradient.png");;</span>
<span class="GAPprompt">gap></span> <span class="GAPinput">F:=ReadImageAsFilteredPureCubicalComplex(file,20);;</span>
<span class="GAPprompt">gap></span> <span class="GAPinput">FF:=ComplementOfFilteredPureCubicalComplex(F);</span>
<span class="GAPprompt">gap></span> <span class="GAPinput">W:=(FiltrationTerm(FF,3));</span>
<span class="GAPprompt">gap></span> <span class="GAPinput">for n in [4..23] do</span>
<span class="GAPprompt">></span> <span class="GAPinput">L:=[];;</span>
<span class="GAPprompt">></span> <span class="GAPinput">for i in [1..PathComponentOfPureComplex(W,0)] do</span>
<span class="GAPprompt">></span> <span class="GAPinput"> P:=PathComponentOfPureComplex(W,i);;</span>
<span class="GAPprompt">></span> <span class="GAPinput"> Q:=ThickeningFiltration(P,150,FiltrationTerm(FF,n));;</span>
<span class="GAPprompt">></span> <span class="GAPinput"> Add(L,Q);;</span>
<span class="GAPprompt">></span> <span class="GAPinput">od;;</span>
<span class="GAPprompt">></span> <span class="GAPinput">W:=Basins(L);</span>
<span class="GAPprompt">></span> <span class="GAPinput">od;</span>
<span class="GAPprompt">gap></span> <span class="GAPinput">C:=PureComplexComplement(W);;</span>
<span class="GAPprompt">gap></span> <span class="GAPinput">T:=PureComplexThickened(C);; C:=ContractedComplex(T,C);; </span>
<span class="GAPprompt">gap></span> <span class="GAPinput">Display(C);</span>
</pre></div>
<p><img src="images/circularGradientSeg.png" align="center" height="250" alt="segmented image"/></p>
<p><a id="X7D2CC9CB85DF1BAF" name="X7D2CC9CB85DF1BAF"></a></p>
<h4>5.8 <span class="Heading">Alternative approaches to computing persistent homology</span></h4>
<p>From any sequence <span class="SimpleMath">X_0 ⊂ X_1 ⊂ X_2 ⊂ ⋯ ⊂ X_T</span> of cellular spaces (such as pure cubical complexes, or cubical complexes, or simplicial complexes, or regular CW complexes) we can construct a filtered chain complex <span class="SimpleMath">C_∗ X_0 ⊂ C_∗ X_1 ⊂ C_∗ X_2 ⊂ ⋯ C_∗ X_T</span>. The induced homology homomorphisms <span class="SimpleMath">H_n(C_∗ X_0, F) → H_n(C_∗ X_1, F) → H_n(C_∗ X_2, F) → ⋯ → H_n(C_∗ X_T, F)</span> with coefficients in a field <span class="SimpleMath">F</span> can be computed by applying an appropriate sequence of elementary row operations to the boundary matrices in the chain complex <span class="SimpleMath">C_∗ X_T⊗ F</span>; the boundary matrices are sparse and are best represented as such; the row operations need to be applied in a fashion that respects the filtration. This method is used in the above examples of persistent homology. The method is not practical when the number of cells in <span class="SimpleMath">X_T</span> is large.</p>
<p>An alternative approach is to construct an admissible discrete vector field on each term <span class="SimpleMath">X_k</span> in the filtration. For each vector field there is a non-regular CW-complex <span class="SimpleMath">Y_k</span> whose cells correspond to the critical cells in <span class="SimpleMath">X_k</span> and for which there is a homotopy equivalence <span class="SimpleMath">X_k≃ Y_k</span>. For each <span class="SimpleMath">k</span> the composite homomorphism <span class="SimpleMath">H_n(C_∗ Y_k, F) stackrel≅→ H_n(C_∗ X_k, F) → H_n(C_∗ X_k+1, F) stackrel≅→ H_n(C_∗ Y_k+1, F)</span> can be computed and the persistent homology can be derived from these homology homomorphisms. This method is implemented in the function <code class="code">PersistentBettiNUmbersAlt(X,n,p)</code> where <span class="SimpleMath">p</span> is the characteristic of the field, <span class="SimpleMath">n</span> is the homology degree, and <span class="SimpleMath">X</span> can be a filtered pure cubical complex, or a filtered simplicial complex, or a filtered regular CW complex, or indeed a filtered chain complex (represented in sparse form). This function incorporates the functions <code class="code">ContractedFilteredPureCubicalComplex(X)</code> and <code class="code">ContractedFilteredRegularComplex(X)</code> which respectively input a filtered pure cubical complex and filtered regular CW-complex and return a filtered complex of the same data type in which each term of the output filtration is a deformation retract of the corresponding term in the input filtration.</p>
<p>In this approach the vector fields on the various spaces <span class="SimpleMath">X_k</span> are completely independent and so the method lends itself to a degree of easy parallelism. This is not incorporated into the current implementation.</p>
<p>As an illustration we consider a synthetic data set <span class="SimpleMath">S</span> consisting of <span class="SimpleMath">3527</span> points sampled, with errors, from an `unknown' manifold <span class="SimpleMath">M</span> in <span class="SimpleMath">R^3</span>. From such a data set one can associate a <span class="SimpleMath">3</span>-dimensional cubical complex <span class="SimpleMath">X_0</span> consisting of one unit cube centred on each (suitably scaled) data point. A visualization of <span class="SimpleMath">X_0</span> is shown below.</p>
<p><img src="images/data.png" align="center" height="300" alt="data cloud"/></p>
<p>Given a pure cubical complex <span class="SimpleMath">X_s</span> we construct <span class="SimpleMath">X_s+1 =X_s ∪ {overline e^3_λ}_λ∈ Λ</span> by adding to <span class="SimpleMath">X_s</span> each closed unit cube <span class="SimpleMath">overline e^3_λ</span> in <span class="SimpleMath">R^3</span> that intersects non-trivially with <span class="SimpleMath">X_s</span>. We construct the filtered cubical complex <span class="SimpleMath">X_∗ ={X_i}_0≤ i≤ 19</span> and compute the persistence matrices <span class="SimpleMath">β_d^∗∗</span> for <span class="SimpleMath">d=0,1,2</span> and for <span class="SimpleMath">Z_2</span> coefficients. The filtered complex <span class="SimpleMath">X_∗</span> is quite large. In particular, the final space <span class="SimpleMath">X_19</span> in the filtration involves <span class="SimpleMath">1092727</span> vertices, <span class="SimpleMath">3246354</span> edges, <span class="SimpleMath">3214836</span> faces of dimension <span class="SimpleMath">2</span> and <span class="SimpleMath">1061208</span> faces of dimension <span class="SimpleMath">3</span>. The usual matrix reduction approach to computing persistent Betti numbers would involve an appropriate row reduction of sparse matrices one of which has over 3 million rows and 3 million columns.</p>
<div class="example"><pre>
<span class="GAPprompt">gap></span> <span class="GAPinput">file:=HapFile("data247.txt");;</span>
<span class="GAPprompt">gap></span> <span class="GAPinput">Read(file);;</span>
<span class="GAPprompt">gap></span> <span class="GAPinput">F:=ThickeningFiltration(T,20);;</span>
<span class="GAPprompt">gap></span> <span class="GAPinput">P:=PersistentBettiNumbersAlt(F,[0,1,2]);;</span>
<span class="GAPprompt">gap></span> <span class="GAPinput">BarCodeCompactDisplay(P);</span>
</pre></div>
<p><img src="images/barcodes123.gif" align="center" height="400" alt="barcodes"/></p>
<p>The barcodes suggest that the data points might have been sampled from a manifold with the homotopy type of a torus.</p>
<p><a id="X86FD0A867EC9E64F" name="X86FD0A867EC9E64F"></a></p>
<h5>5.8-1 <span class="Heading">Non-trivial cup product</span></h5>
<p>Of course, a wedge <span class="SimpleMath">S^2∨ S^1∨ S^1</span> has the same homology as the torus <span class="SimpleMath">S^1× S^1</span>. By establishing that a 'noise free' model for our data points, say the 10-th term <span class="SimpleMath">X_10</span> in the filtration, has a non-trivial cup product <span class="SimpleMath">∪: H^1(X_10, Z) × H^1(X_10, Z) → H^2(X_10, Z)</span> we can eliminate <span class="SimpleMath">S^2∨ S^1∨ S^1</span> as a candidate from which the data was sampled.</p>
<div class="example"><pre>
<span class="GAPprompt">gap></span> <span class="GAPinput">X10:=RegularCWComplex(FiltrationTerm(F,10));;</span>
<span class="GAPprompt">gap></span> <span class="GAPinput">cup:=LowDimensionalCupProduct(X10);;</span>
<span class="GAPprompt">gap></span> <span class="GAPinput">cup([1,0],[0,1]);</span>
[ 1 ]
</pre></div>
<p><a id="X783EF0F17B629C46" name="X783EF0F17B629C46"></a></p>
<h5>5.8-2 <span class="Heading">Explicit homology generators</span></h5>
<p>It could be desirable to obtain explicit representatives of the persistent homology generators that "persist" through a significant sequence of filtration terms. There are two such generators in degree <span class="SimpleMath">1</span> and one such generator in degree <span class="SimpleMath">2</span>. The explicit representatives in degree <span class="SimpleMath">n</span> could consist of an inclusion of pure cubical complexes <span class="SimpleMath">Y_n ⊂ X_10</span> for which the incuced homology homomorphism <span class="SimpleMath">H_n(Y_n, Z) → H_n(X_10, Z)</span> is an isomorphism, and for which <span class="SimpleMath">Y_n</span> is minimal in the sense that its homotopy type changes if any one or more of its top dimensional cells are removed. Ideally the space <span class="SimpleMath">Y_n</span> should be "close to the original dataset" <span class="SimpleMath">X_0</span>. The following commands first construct an explicit degree <span class="SimpleMath">2</span> homology generator representative <span class="SimpleMath">Y_2⊂ X_10</span> where <span class="SimpleMath">Y_2</span> is homotopy equivalent to <span class="SimpleMath">X_10</span>. They then construct an explicit degree <span class="SimpleMath">1</span> homology generators representative <span class="SimpleMath">Y_1⊂ X_10</span> where <span class="SimpleMath">Y_1</span> is homotopy equivalent to a wedge of two circles. The final command displays the homology generators representative <span class="SimpleMath">Y_1</span>.</p>
<div class="example"><pre>
<span class="GAPprompt">gap></span> <span class="GAPinput">Y2:=FiltrationTerm(F,10);; </span>
<span class="GAPprompt">gap></span> <span class="GAPinput">for t in Reversed([1..9]) do</span>
<span class="GAPprompt">></span> <span class="GAPinput">Y2:=ContractedComplex(Y2,FiltrationTerm(F,t));</span>
<span class="GAPprompt">></span> <span class="GAPinput">od;</span>
<span class="GAPprompt">gap></span> <span class="GAPinput">Y2:=ContractedComplex(Y2);;</span>
<span class="GAPprompt">gap></span> <span class="GAPinput">Size(FiltrationTerm(F,10));</span>
918881
<span class="GAPprompt">gap></span> <span class="GAPinput">Size(Y2); </span>
61618
<span class="GAPprompt">gap></span> <span class="GAPinput">Y1:=PureComplexDifference(Y2,PureComplexRandomCell(Y2));;</span>
<span class="GAPprompt">gap></span> <span class="GAPinput">Y1:=ContractedComplex(Y1);;</span>
<span class="GAPprompt">gap></span> <span class="GAPinput">Size(Y1);</span>
474
<span class="GAPprompt">gap></span> <span class="GAPinput">Display(Y1);</span>
</pre></div>
<p><img src="images/cubicaltorusgens.png" align="center" height="200" alt="first homology generators"/></p>
<p><a id="X80D0D8EB7BCD05E9" name="X80D0D8EB7BCD05E9"></a></p>
<h4>5.9 <span class="Heading">Knotted proteins</span></h4>
<p>The <span class="URL"><a href="https://www.rcsb.org/">Protein Data Bank</a></span> contains a wealth of data which can be investigated with respect to knottedness. Information on a particular protein can be downloaded as a .pdb file. Each protein consists of one or more chains of amino acids and the file gives 3-dimensional Euclidean coordinates of the atoms in amino acids. Each amino acid has a unique "alpha carbon" atom (labelled as "CA" in the pdb file). A simple 3-dimensional curve, the <em>protein backbone</em>, can be constructed through the sequence of alpha carbon atoms. Typically the ends of the protein backbone lie near the "surface" of the protein and can be joined by a path outside of the protein to obtain a simple closed curve in Euclidean 3-space.</p>
<p>The following command reads in the pdb file for the T.thermophilus 1V2X protein, which consists of a single chain of amino acids, and uses Asymptote software to produce an interactive visualization of its backbone. A path joining the end vertices of the backbone is displayed in blue.</p>
<div class="example"><pre>
<span class="GAPprompt">gap></span> <span class="GAPinput">file:=HapFile("data1V2X.pdb");;</span>
<span class="GAPprompt">gap></span> <span class="GAPinput">DisplayPDBfile(file);</span>
</pre></div>
<p><img src="images/1v2x.gif" align="center" height="500" alt="a protein backbone"/></p>
<p>The next command reads in the pdb file for the T.thermophilus 1V2X protein and represents it as a <span class="SimpleMath">3</span>-dimensional pure cubical complex <span class="SimpleMath">K</span>. A resolution of <span class="SimpleMath">r=5</span> is chosen and this results in a representation as a subcomplex <span class="SimpleMath">K</span> of an ambient rectangular box of volume equal to <span class="SimpleMath">184× 186× 294</span> unit cubes. The complex <span class="SimpleMath">K</span> should have the homotopy type of a circle and the protein backbone is a 1-dimenional curve that should lie in <span class="SimpleMath">K</span>. The final command displays <span class="SimpleMath">K</span>.</p>
<div class="example"><pre>
<span class="GAPprompt">gap></span> <span class="GAPinput">r:=5;;</span>
<span class="GAPprompt">gap></span> <span class="GAPinput">K:=ReadPDBfileAsPureCubicalComplex(file,r);; </span>
<span class="GAPprompt">gap></span> <span class="GAPinput">K:=ContractedComplex(K);;</span>
<span class="GAPprompt">gap></span> <span class="GAPinput">K!.properties;</span>
[ [ "dimension", 3 ], [ "arraySize", [ 184, 186, 294 ] ] ]
<span class="GAPprompt">gap></span> <span class="GAPinput">Display(K);</span>
</pre></div>
<p><img src="images/1v2xcubical.gif" align="center" height="500" alt="pure cubical complex representing a protein backbone"/></p>
<p>Next we create a filtered pure cubical complex by repeatedly thickening <span class="SimpleMath">K</span>. We perform <span class="SimpleMath">15</span> thickenings, each thickening being a term in the filtration. The <span class="SimpleMath">β_1</span> barcode for the filtration is displayed. This barcode is a descriptor for the geometry of the protein. For current purposes it suffices to note that the first few terms of the filtration have first homology equal to that of a circle. This indicates that the Euclidean coordinates in the pdb file robustly determine some knot.</p>
<div class="example"><pre>
<span class="GAPprompt">gap></span> <span class="GAPinput">F:=ThickeningFiltration(K,15);;</span>
<span class="GAPprompt">gap></span> <span class="GAPinput">F:=FilteredPureCubicalComplexToCubicalComplex(F);;</span>
<span class="GAPprompt">gap></span> <span class="GAPinput">F:=FilteredCubicalComplexToFilteredRegularCWComplex(F);;</span>
<span class="GAPprompt">gap></span> <span class="GAPinput">P:=PersistentBettiNumbersAlt(F,1);;</span>
<span class="GAPprompt">gap></span> <span class="GAPinput">BarCodeCompactDisplay(P);</span>
</pre></div>
<p><img src="images/1v2xbarcode.gif" align="center" height="500" alt="barcode"/></p>
<p>The next commands compute a presentation for the fundamental group <span class="SimpleMath">π_1( R^3∖ K)</span> and the Alexander polynomial for the knot. This is the same Alexander polynomial as for the trefoil knot. Also, Tietze transformations can be used to see that the fundamental group is the same as for the trefoil knot.</p>
<div class="example"><pre>
<span class="GAPprompt">gap></span> <span class="GAPinput">C:=PureComplexComplement(K);;</span>
<span class="GAPprompt">gap></span> <span class="GAPinput">C:=ContractedComplex(C);;</span>
<span class="GAPprompt">gap></span> <span class="GAPinput">G:=FundamentalGroup(C);;</span>
<span class="GAPprompt">gap></span> <span class="GAPinput">GeneratorsOfGroup(G);</span>
[ f1, f2 ]
<span class="GAPprompt">gap></span> <span class="GAPinput">RelatorsOfFpGroup(G);</span>
[ f2*f1^-1*f2^-1*f1^-1*f2*f1 ]
<span class="GAPprompt">gap></span> <span class="GAPinput">AlexanderPolynomial(G);</span>
x_1^2-x_1+1
</pre></div>
<p><a id="X87AF06677F05C624" name="X87AF06677F05C624"></a></p>
<h4>5.10 <span class="Heading">Random simplicial complexes</span></h4>
<p>For a positive integer <span class="SimpleMath">n</span> and probability <span class="SimpleMath">p</span> we denote by <span class="SimpleMath">Y(n,p)</span> the <em>Linial-Meshulam random simplicial 2-complex</em>. Its <span class="SimpleMath">1</span>-skeleton is the complete graph on <span class="SimpleMath">n</span> vertices; each possible <span class="SimpleMath">2</span>-simplex is included independently with probability <span class="SimpleMath">p</span>.</p>
<p>The following commands first compute the number <span class="SimpleMath">h_i</span> of non-trivial cyclic summands in <span class="SimpleMath">H_i(Y(100,p), Z)</span> for a range of probabilities <span class="SimpleMath">p</span> and <span class="SimpleMath">i=1,2</span> and then produce a plot of <span class="SimpleMath">h_i</span> versus <span class="SimpleMath">p</span>. The plot for <span class="SimpleMath">h_1</span> is red and the plot for <span class="SimpleMath">h_2</span> is blue. A plot for the Euler characteristic <span class="SimpleMath">1-h_1+h_2</span> is shown in green.</p>
<div class="example"><pre>
<span class="GAPprompt">gap></span> <span class="GAPinput">L:=[];;M:=[];;</span>
<span class="GAPprompt">gap></span> <span class="GAPinput">for p in [1..100] do</span>
<span class="GAPprompt">></span> <span class="GAPinput">K:=RegularCWComplex(RandomSimplicialTwoComplex(100,p/1000));;</span>
<span class="GAPprompt">></span> <span class="GAPinput">h1:=Length(Homology(K,1));;</span>
<span class="GAPprompt">></span> <span class="GAPinput">h2:=Length(Homology(K,2));;</span>
<span class="GAPprompt">></span> <span class="GAPinput">Add(L, [1.0*(p/1000),h1,"red"]);</span>
<span class="GAPprompt">></span> <span class="GAPinput">Add(L, [1.0*(p/1000),h2,"blue"]);</span>
<span class="GAPprompt">></span> <span class="GAPinput">Add(M, [1.0*(p/1000),1-h1+h2,"green"]);</span>
<span class="GAPprompt">></span> <span class="GAPinput">od;</span>
<span class="GAPprompt">gap></span> <span class="GAPinput">ScatterPlot(L);</span>
<span class="GAPprompt">gap></span> <span class="GAPinput">ScatterPlot(M);</span>
</pre></div>
<p><img src="tutex/graph4.7.png" align="center" height="500" alt="a graph"/> <img src="tutex/graph4.77.png" align="center" height="500" alt="a graph"/></p>
<p>From this plot it seems that there is a <em>phase change threshold</em> at around <span class="SimpleMath">p=0.025</span>. An inspection of the first homology groups <span class="SimpleMath">H_1(Y(100,p), Z)</span> shows that in most cases there is no torsion. However, around the threshold some of the complexes do have torsion in their first homology.</p>
<p>Similar commands for <span class="SimpleMath">Y(75,p)</span> suggest a phase transition at around <span class="SimpleMath">p=0.035</span> in this case. The following commands compute <span class="SimpleMath">H_1(Y(75,p), Z)</span> for <span class="SimpleMath">900</span> random <span class="SimpleMath">2</span>-complexes with <span class="SimpleMath">p</span> in a small interval around <span class="SimpleMath">0.035</span> and, in each case where there is torsion, the torsion coefficients are stored in a list. The final command prints these lists -- all but one of which are of length <span class="SimpleMath">1</span>. For example, there is one <span class="SimpleMath">2</span>-dimensional simplicial complex on <span class="SimpleMath">75</span> vertices whose first homology contains the summand <span class="SimpleMath">Z_107879661870516800665161182578823128</span>. The largest prime factor is <span class="SimpleMath">80555235907994145009690263</span> occuring in the summand <span class="SimpleMath">Z_259837760616287294231081766978855</span>.</p>
<div class="example"><pre>
<span class="GAPprompt">gap></span> <span class="GAPinput">torsion:=function(n,p)</span>
<span class="GAPprompt">></span> <span class="GAPinput">local H, Y;</span>
<span class="GAPprompt">></span> <span class="GAPinput">Y:=RegularCWComplex(RandomSimplicialTwoComplex(n,p));</span>
<span class="GAPprompt">></span> <span class="GAPinput">H:=Homology(Y,1);</span>
<span class="GAPprompt">></span> <span class="GAPinput">H:=Filtered(H,x->not x=0);</span>
<span class="GAPprompt">></span> <span class="GAPinput">return H;</span>
<span class="GAPprompt">></span> <span class="GAPinput">end;</span>
function( n, p ) ... end
<span class="GAPprompt">gap></span> <span class="GAPinput">L:=[];;for n in [73000..73900] do</span>
<span class="GAPprompt">></span> <span class="GAPinput">t:=torsion(75,n/2000000); </span>
<span class="GAPprompt">></span> <span class="GAPinput">if not t=[] then Add(L,t); fi;</span>
<span class="GAPprompt">></span> <span class="GAPinput">od;</span>
<span class="GAPprompt">gap></span> <span class="GAPinput">Display(L);</span>
[ [ 2 ],
[ 26 ],
[ 259837760616287294231081766978855 ],
[ 2 ],
[ 3 ],
[ 2 ],
[ 2761642698060127444812143568 ],
[ 2626355281010974663776273381976 ],
[ 2 ],
[ 3 ],
[ 33112382751264894819430785350 ],
[ 16 ],
[ 4 ],
[ 3 ],
[ 2 ],
[ 3 ],
[ 2 ],
[ 85234949999183888967763100590977 ],
[ 2 ],
[ 24644196130785821107897718662022 ],
[ 2, 2 ],
[ 2 ],
[ 416641662889025645492982468 ],
[ 41582773001875039168786970816 ],
[ 2 ],
[ 75889883165411088431747730 ],
[ 33523474091636554792305315165 ],
[ 107879661870516800665161182578823128 ],
[ 5588265814409119568341729980 ],
[ 2 ],
[ 5001457249224115878015053458 ],
[ 10 ],
[ 12 ],
[ 2 ],
[ 2 ],
[ 3 ],
[ 7757870243425246987971789322 ],
[ 8164648856993269673396613497412 ],
[ 2 ] ]
</pre></div>
<p><a id="X875EE92F7DBA1E27" name="X875EE92F7DBA1E27"></a></p>
<h4>5.11 <span class="Heading">Computing homology of a clique complex (Vietoris-Rips complex) </span></h4>
<p>Topological data analysis provides one motivation for wanting to compute the homology of a clique complex. Consider for instance the cloud of data points shown in Example <a href="chap5.html#X7D512DA37F789B4C"><span class="RefLink">5.2-1</span></a>. This data is a set <span class="SimpleMath">S</span> of 400 points in the plane. Let <span class="SimpleMath">Γ</span> be the graph with vertex set <span class="SimpleMath">S</span> and with two vertices joined by an edge if they lie within a Euclidean distance of 40 of each other. The clique complex <span class="SimpleMath">K=K(Γ)</span> could be studied to see what it reveals about the data. The following commands construct <span class="SimpleMath">K</span> and show that it is a 23-dimensional simplicial complex consisting of a total of 36191976 simplices.</p>
<div class="example"><pre>
<span class="GAPprompt">gap></span> <span class="GAPinput">file:=HapFile("data134.txt");; </span>
<span class="GAPprompt">gap></span> <span class="GAPinput">Read(file);</span>
<span class="GAPprompt">gap></span> <span class="GAPinput">A:=VectorsToSymmetricMatrix(S,EuclideanApproximatedMetric);;</span>
<span class="GAPprompt">gap></span> <span class="GAPinput">threshold:=40;; </span>
<span class="GAPprompt">gap></span> <span class="GAPinput">grph:=SymmetricMatrixToGraph(A,threshold);;</span>
<span class="GAPprompt">gap></span> <span class="GAPinput">dimension_cap:=100;; </span>
<span class="GAPprompt">gap></span> <span class="GAPinput">K:=CliqueComplex(grph,dimension_cap);</span>
Simplicial complex of dimension 23.
<span class="GAPprompt">gap></span> <span class="GAPinput">Size(K);</span>
36191976
</pre></div>
<p>The computation of the homology of this clique complex <span class="SimpleMath">K</span> is a challenge because of its size. If we are only interested in <span class="SimpleMath">K</span> up to homotopy then we could try to modify the graph <span class="SimpleMath">Γ</span> in such a way that the homotopy type of the clique complex is unchanged but the size of the clique complex is reduced. This is done in the following commands, producing a smaller <span class="SimpleMath">19</span>-dimensional simplicial complex <span class="SimpleMath">K</span> with 4180652 simplices.</p>
<div class="example"><pre>
<span class="GAPprompt">gap></span> <span class="GAPinput">ContractGraph(grph);;</span>
<span class="GAPprompt">gap></span> <span class="GAPinput">dimension_cap:=100;; </span>
<span class="GAPprompt">gap></span> <span class="GAPinput">K:=CliqueComplex(grph,dimension_cap);</span>
Simplicial complex of dimension 19.
<span class="GAPprompt">gap></span> <span class="GAPinput">Size(K);</span>
4180652
</pre></div>
<p>To compute the homology of <span class="SimpleMath">K</span> in degrees <span class="SimpleMath">0</span> to <span class="SimpleMath">5</span> say, we could represent <span class="SimpleMath">K</span> as a regular CW-complex <span class="SimpleMath">Y</span> and then compute the homology of <span class="SimpleMath">Y</span> as follows. The homology <span class="SimpleMath">H_n(K)= Z</span> for <span class="SimpleMath">n=0,1</span> and <span class="SimpleMath">H_n(K)= 0</span> for <span class="SimpleMath">n=2,3,4,5</span> is consistent with the data having been sampled from a space with the homotopy type of a circle.</p>
<div class="example"><pre>
<span class="GAPprompt">gap></span> <span class="GAPinput">Y:=RegularCWComplex(K);</span>
Regular CW-complex of dimension 19
<span class="GAPprompt">gap></span> <span class="GAPinput">Homology(Y,0);</span>
[ 0 ]
<span class="GAPprompt">gap></span> <span class="GAPinput">Homology(Y,1);</span>
[ 0 ]
<span class="GAPprompt">gap></span> <span class="GAPinput">Homology(Y,2);</span>
[ ]
<span class="GAPprompt">gap></span> <span class="GAPinput">Homology(Y,3);</span>
[ ]
<span class="GAPprompt">gap></span> <span class="GAPinput">Homology(Y,4);</span>
[ ]
<span class="GAPprompt">gap></span> <span class="GAPinput">Homology(Y,5)</span>
[ ]
</pre></div>
<div class="chlinkprevnextbot"> <a href="chap0.html">[Top of Book]</a> <a href="chap0.html#contents">[Contents]</a> <a href="chap4.html">[Previous Chapter]</a> <a href="chap6.html">[Next Chapter]</a> </div>
<div class="chlinkbot"><span class="chlink1">Goto Chapter: </span><a href="chap0.html">Top</a> <a href="chap1.html">1</a> <a href="chap2.html">2</a> <a href="chap3.html">3</a> <a href="chap4.html">4</a> <a href="chap5.html">5</a> <a href="chap6.html">6</a> <a href="chap7.html">7</a> <a href="chap8.html">8</a> <a href="chap9.html">9</a> <a href="chap10.html">10</a> <a href="chap11.html">11</a> <a href="chap12.html">12</a> <a href="chap13.html">13</a> <a href="chap14.html">14</a> <a href="chap15.html">15</a> <a href="chap16.html">16</a> <a href="chapBib.html">Bib</a> <a href="chapInd.html">Ind</a> </div>
<hr />
<p class="foot">generated by <a href="https://www.math.rwth-aachen.de/~Frank.Luebeck/GAPDoc">GAPDoc2HTML</a></p>
</body>
</html>
|
