David Dumas - Python for mathematical visualization: a four-dimensional case study

This is a talk about creating pictures of a mathematical object---specifically, a 4-dimensional fractal "dust" that has been the subject of mathematical research in hyperbolic geometry since the 1980s. In the end this is accomplished using a little algebra, a little geometry, and a healthy dose of Python.

That is, I will present a case study of using Python in several aspects of a mathematical visualization project, from the computation itself, to transforming and converting data, and finally for scripting the process of generating the images. Along the way I'll explain how Python's convenient idioms and containers (e.g. sets and set comprehensions) are a good fit for some of the algebraic and geometric questions that come up, how Scipy and Numpy enable fast numerical calculations, and how Python's strength as a language for scripting and automation allows easy orchestration of rendering of still images and frames of animations.

The mathematical visualization project we describe is a collaboration with François Guéritaud (Université de Lille).

https://us.pycon.org/2017/schedule/presentation/746/

Transcript

1. Python for Mathematical Visualization A Four-Dimensional Case Study David Dumas University of Illinois at Chicago May 20, 2017
2. Goal Describe a mathematical visualization project that uses Python. Highlight characteristics of Python that make it suitable for this application.
3. The PML Visualization Project Joint work with François Guéritaud http://dumas.io/PML
4. Steps dumas.io/PML Enumerate the objects Calculate coordinates Render an image
5. Topological enumeration dumas.io/PML START FINISH
6. Topological enumeration dumas.io/PML START FINISH Obstacle
7. Topological enumeration dumas.io/PML START FINISH
8. Topological enumeration dumas.io/PML START FINISH
9. Topological enumeration dumas.io/PML START FINISH
10. Topological enumeration dumas.io/PML START FINISH
11. Topological enumeration dumas.io/PML START FINISH
12. Topological enumeration dumas.io/PML START FINISH
13. Topological enumeration dumas.io/PML START FINISH Not a simple path
14. Topological enumeration dumas.io/PML START FINISH
15. Topological enumeration dumas.io/PML START FINISH Simple path space
16. Topological enumeration dumas.io/PML Simple loops on the five-holed sphere Holes
17. Topological enumeration dumas.io/PML Simple loops on the five-holed sphere START & FINISH
18. Topological enumeration dumas.io/PML Simple loops on the five-holed sphere START & FINISH
19. Topological enumeration dumas.io/PML Simple loops on the five-holed sphere
20. Topological enumeration dumas.io/PML Simple loops on the five-holed sphere
21. Topological enumeration dumas.io/PML Simple loops on the five-holed sphere
22. Topological enumeration dumas.io/PML Simple loops on the five-holed sphere
23. Topological enumeration dumas.io/PML Simple loops on the five-holed sphere
24. Topological enumeration dumas.io/PML Simple loops on the five-holed sphere Holes
25. Topological enumeration dumas.io/PML Simple loops on the five-holed sphere
26. Topological enumeration dumas.io/PML Simple loops on the five-holed sphere
27. Topological enumeration dumas.io/PML Simple loops on the five-holed sphere 2 3
28. Topological enumeration dumas.io/PML Nonperipheral simple closed curves on the five-holed sphere 2 3
29. Topological enumeration dumas.io/PML Nonperipheral simple closed curves on the five-holed sphere (NSCC) 2 3
30. Topological enumeration dumas.io/PML Nonperipheral simple closed curves on the five-holed sphere (NSCC)
31. Topological enumeration dumas.io/PML Nonperipheral simple closed curves on the five-holed sphere (NSCC)
32. Topological enumeration dumas.io/PML Nonperipheral simple closed curves on the five-holed sphere (NSCC)
33. Topological enumeration dumas.io/PML Nonperipheral simple closed curves on the five-holed sphere (NSCC)
34. Topological enumeration dumas.io/PML Nonperipheral simple closed curves on the five-holed sphere (NSCC)
35. Topological enumeration dumas.io/PML Nonperipheral simple closed curves on the five-holed sphere (NSCC)
36. Topological enumeration dumas.io/PML Nonperipheral simple closed curves on the five-holed sphere (NSCC)
37. Topological enumeration dumas.io/PML Nonperipheral simple closed curves on the five-holed sphere (NSCC)
38. Topological enumeration dumas.io/PML Nonperipheral simple closed curves on the five-holed sphere (NSCC)
39. Generators dumas.io/PML a
40. Generators dumas.io/PML a b
41. Generators dumas.io/PML a b c
42. Generators dumas.io/PML a b c d
43. Generators dumas.io/PML a b c d e
44. Generators dumas.io/PML A B C D E
45. Generators dumas.io/PML a b c d e
46. Encoding curves in strings dumas.io/PML ab a b c d e
47. Encoding curves in strings dumas.io/PML ab a b c d e
48. Encoding curves in strings dumas.io/PML ab a b c d e
49. Encoding curves in strings dumas.io/PML ab a b c d e
50. Encoding curves in strings dumas.io/PML ab a b c d e
51. Encoding curves in strings dumas.io/PML bc a b c d e
52. Encoding curves in strings dumas.io/PML ac a b c d e
53. Encoding curves in strings dumas.io/PML CA a b c d e
54. Encoding curves in strings dumas.io/PML aC a b c d e
55. Encoding curves in strings dumas.io/PML aC a b c d e
56. Encoding curves in strings dumas.io/PML abDcd a b c d e
57. Naïve generator dumas.io/PML All words in alphabet 'abcdeABCDE'. from itertools import product N = 5 # Word length for w in product('abcdeABCDE',repeat=N): print(''.join(w))
58. Naïve generator dumas.io/PML Simple words in alphabet 'abcdeABCDE'. from itertools import product N = 5 # Word length for w in product('abcdeABCDE',repeat=N): if is_simple_curve(w): print(''.join(w))
59. Naïve generator dumas.io/PML Simple words in alphabet 'abcdeABCDE'. from itertools import product N = 5 # Word length for w in product('abcdeABCDE',repeat=N): if is_simple_curve(w): print(''.join(w)) # TODO: Implement is_simple_curve()
60. Exponential haystack dumas.io/PML aa, ab, ac, ad, ae, aB, aC, aD, aE, ba, bb, bc, bd, be, bA, bC, bD, bE, ca, cb, cc, cd, ce, cA, cB, cD, cE, da, db, dc, dd, de, dA, dB, dC, dE, ea, eb, ec, ed, ee, eA, eB, eC, eD, Ab, Ac, Ad, Ae, AA, AB, AC, AD, AE, Ba, Bc, Bd, Be, BA, BB, BC, BD, BE, Ca, Cb, Cd, Ce, CA, CB, CC, CD, CE, Da, Db, Dc, De, DA, DB, DC, DD, DE, Ea, Eb, Ec, Ed, EA, EB, EC, ED, EE
61. Exponential haystack dumas.io/PML aaa, aab, aac, aad, aae, aaB, aaC, aaD, aaE, aba, abb, abc, abd, abe, abA, abC, abD, abE, aca, acb, acc, acd, ace, acA, acB, acD, acE, ada, adb, adc, add, ade, adA, adB, adC, adE, aea, aeb, aec, aed, aee, aeA, aeB, aeC, aeD, aBa, aBc, aBd, aBe, aBA, aBB, aBC, aBD, aBE, aCa, aCb, aCd, aCe, aCA, aCB, aCC, aCD, aCE, aDa, aDb, aDc, aDe, aDA, aDB, aDC, aDD, aDE, aEa, aEb, aEc, aEd, aEA, aEB, aEC, aED, aEE, baa, bab, bac, bad, bae, baB, baC, baD, baE, bba, bbb, bbc, bbd, bbe, bbA, bbC, bbD, bbE, bca, bcb, bcc, bcd, bce, bcA, bcB, bcD, bcE, bda, bdb, bdc, bdd, bde, bdA, bdB, bdC, bdE, bea, beb, bec, bed, bee, beA, beB, beC, beD, bAb, bAc, bAd, bAe, bAA, bAB, bAC, bAD, bAE, bCa, bCb, bCd, bCe, bCA, bCB, bCC, bCD, bCE, bDa, bDb, bDc, bDe, bDA, bDB, bDC, bDD, bDE, bEa, bEb, bEc, bEd, bEA, bEB, bEC, bED, bEE, caa, cab, cac, cad, cae, caB, caC, caD, caE, cba, cbb, cbc, cbd, cbe, cbA, cbC, cbD, cbE, cca, ccb, ccc, ccd, cce, ccA, ccB, ccD, ccE, cda, cdb, cdc, cdd, cde, cdA, cdB, cdC, cdE, cea, ceb, cec, ced, cee, ceA, ceB, ceC, ceD, cAb, cAc, cAd, cAe, cAA, cAB, cAC, cAD, cAE, cBa, cBc, cBd, cBe, cBA, cBB, cBC, cBD, cBE, cDa, cDb, cDc, cDe, cDA, cDB, cDC, cDD, cDE, cEa, cEb, cEc, cEd, cEA, cEB, cEC, cED, cEE, daa, dab, dac, dad, dae, daB, daC, daD, daE, dba, dbb, dbc, dbd, dbe, dbA, dbC, dbD, dbE, dca, dcb, dcc, dcd, dce, dcA, dcB, dcD, dcE, dda, ddb, ddc, ddd, dde, ddA, ddB, ddC, ddE, dea, deb, dec, ded, dee, deA, deB, deC, deD, dAb, dAc, dAd, dAe, dAA, dAB, dAC, dAD, dAE, dBa, dBc, dBd, dBe, dBA, dBB, dBC, dBD, dBE, dCa, dCb, dCd, dCe, dCA, dCB, dCC, dCD, dCE, dEa, dEb, dEc, dEd, dEA, dEB, dEC, dED, dEE, eaa, eab, eac, ead, eae, eaB, eaC, eaD, eaE, eba, ebb, ebc, ebd, ebe, ebA, ebC, ebD, ebE, eca, ecb, ecc, ecd, ece, ecA, ecB, ecD, ecE, eda, edb, edc, edd, ede, edA, edB, edC, edE, eea, eeb, eec, eed, eee, eeA, eeB, eeC, eeD, eAb, eAc, eAd, eAe, eAA, eAB, eAC, eAD, eAE, eBa, eBc, eBd, eBe, eBA, eBB, eBC, eBD, eBE, eCa, eCb, eCd, eCe, eCA, eCB, eCC, eCD, eCE, eDa, eDb, eDc, eDe, eDA, eDB, eDC, eDD, eDE, Aba, Abb, Abc, Abd, Abe, AbA, AbC, AbD, AbE, Aca, Acb, Acc, Acd, Ace, AcA, AcB, AcD, AcE, Ada, Adb, Adc, Add, Ade, AdA, AdB, AdC, AdE, Aea, Aeb, Aec, Aed, Aee, AeA, AeB, AeC, AeD, AAb, AAc, AAd, AAe, AAA, AAB, AAC, AAD, AAE, ABa, ABc, ABd, ABe, ABA, ABB, ABC, ABD, ABE, ACa, ACb, ACd, ACe, ACA, ACB, ACC, ACD, ACE, ADa, ADb, ADc, ADe, ADA, ADB, ADC, ADD, ADE, AEa, AEb, AEc, AEd, AEA, AEB, AEC, AED, AEE, Baa, Bab, Bac, Bad, Bae, BaB, BaC, BaD, BaE, Bca, Bcb, Bcc, Bcd, Bce, BcA, BcB, BcD, BcE, Bda, Bdb, Bdc, Bdd, Bde, BdA, BdB, BdC, BdE, Bea, Beb, Bec, Bed, Bee, BeA, BeB, BeC, BeD, BAb, BAc, BAd, BAe, BAA, BAB, BAC, BAD, BAE, BBa, BBc, BBd, BBe, BBA, BBB, BBC, BBD, BBE, BCa, BCb, BCd, BCe, BCA, BCB, BCC, BCD, BCE, BDa, BDb, BDc, BDe, BDA, BDB, BDC, BDD, BDE, BEa, BEb, BEc, BEd, BEA, BEB, BEC, BED, BEE, Caa, Cab, Cac, Cad, Cae, CaB, CaC, CaD, CaE, Cba, Cbb, Cbc, Cbd, Cbe, CbA, CbC, CbD, CbE, Cda, Cdb, Cdc, Cdd, Cde, CdA, CdB, CdC, CdE, Cea, Ceb, Cec, Ced, Cee, CeA, CeB, CeC, CeD, CAb, CAc, CAd, CAe, CAA, CAB, CAC, CAD, CAE, CBa, CBc, CBd, CBe, CBA, CBB, CBC, CBD, CBE, CCa, CCb, CCd, CCe, CCA, CCB, CCC, CCD, CCE, CDa, CDb, CDc, CDe, CDA, CDB, CDC, CDD, CDE, CEa, CEb, CEc, CEd, CEA, CEB, CEC, CED, CEE, Daa, Dab, Dac, Dad, Dae, DaB, DaC, DaD, DaE, Dba, Dbb, Dbc, Dbd, Dbe, DbA, DbC, DbD, DbE, Dca, Dcb, Dcc, Dcd, Dce, DcA, DcB, DcD, DcE, Dea, Deb, Dec, Ded, Dee, DeA, DeB, DeC, DeD, DAb, DAc, DAd, DAe, DAA, DAB, DAC, DAD, DAE, DBa, DBc, DBd, DBe, DBA, DBB, DBC, DBD, DBE, DCa, DCb, DCd, DCe, DCA, DCB, DCC, DCD, DCE, DDa, DDb, DDc, DDe, DDA, DDB, DDC, DDD, DDE, DEa, DEb, DEc, DEd, DEA, DEB, DEC, DED, DEE, Eaa, Eab, Eac, Ead, Eae, EaB, EaC, EaD, EaE, Eba, Ebb, Ebc, Ebd, Ebe, EbA, EbC, EbD, EbE, Eca, Ecb, Ecc, Ecd, Ece, EcA, EcB, EcD, EcE, Eda, Edb, Edc, Edd, Ede, EdA, EdB, EdC, EdE, EAb, EAc, EAd, EAe, EAA, EAB, EAC, EAD, EAE, EBa, EBc, EBd, EBe, EBA, EBB, EBC, EBD, EBE, ECa, ECb, ECd, ECe, ECA, ECB, ECC, ECD, ECE, EDa, EDb, EDc, EDe, EDA, EDB, EDC, EDD, EDE, EEa, EEb, EEc, EEd, EEA, EEB, EEC, EED, EEE
62. Exponential haystack dumas.io/PML
63. Twisting dumas.io/PML bc a b c d e cd
64. Twisting dumas.io/PML bc a b c d e cd
65. Twisting dumas.io/PML bc a b c d e ? cd
66. Twisting dumas.io/PML bc a b c d e ? cd
67. Twisting dumas.io/PML bc a b c d e ? cd
68. Twisting dumas.io/PML bc a b c d e bd Tcd
69. Twisting dumas.io/PML bc a b c d e bd Tcd bDcd 2
70. Twisting is search & replace dumas.io/PML Tcd is equivalent to substituting: a → a b → b c → d d → Dcd e → e
71. Twisting is search & replace dumas.io/PML Tab Tbc Tcd Tde Tea a → b a a a Aea b → Bab c b b b c → c Cbc d c c d → d d Dcd e d e → e e e Ede a
72. Twists generate everything dumas.io/PML Fact: Any nonperipheral simple closed curve can be transformed to any other one by applying twists Tab , Tbc , Tcd , Tde , Tea and their inverses.
73. Twist-based generator dumas.io/PML seed curves
74. Twist-based generator dumas.io/PML twists seed curves
75. Twist-based generator dumas.io/PML twists frontier
76. Twist-based generator dumas.io/PML twists frontier
77. Twist-based generator dumas.io/PML twists frontier
78. Twist-based generator dumas.io/PML curves
79. Twist-based generator dumas.io/PML
80. Twist-based generator dumas.io/PML
81. Twist-based generator dumas.io/PML
82. Twist-based generator dumas.io/PML
83. Twist-based generator dumas.io/PML
84. Twist-based generator dumas.io/PML
85. Twist-based generator dumas.io/PML curves
86. Twist-based generator dumas.io/PML Python’s set type stores collections of distinct elements and supports efficient boolean operations. depth = 5 twists = { Tab, Tab_inv, Tbc, Tbc_inv } # etc curves = set() frontier = {'ab', 'bc', 'cd', 'de', 'ea'} for _ in range(depth): latest = \ { T(x) for x in frontier for T in twists } frontier = latest.difference(curves) curves.update(frontier)
87. Coordinates dumas.io/PML How to visualize the entire set of NSCCs? Assign coordinates.
88. Coordinates dumas.io/PML How to visualize the entire set of NSCCs? Assign coordinates. In 1986, W. P. Thurston described a way to associate a 4-dimensional vector to each NSCC on the five-holed sphere.
89. Coordinates dumas.io/PML How to visualize the entire set of NSCCs? Assign coordinates. In 1986, W. P. Thurston described a way to associate a 4-dimensional vector to each NSCC on the five-holed sphere. The resulting cloud of points densely fills a 3-dimensional hypersurface in 4-space. This hypersurface is PML.
90. Length & dlength dumas.io/PML Curve x Tiling matrix (x) length L(x)
91. Length & dlength dumas.io/PML Curve x Tiling matrix (x) length L(x) 6.467 . . .
92. Length & dlength dumas.io/PML Curve x Tiling matrix (x) length L(x) 6.467 . . . 1 1 1 1
93. Length & dlength dumas.io/PML Curve x Tiling matrix (x) length L(x) 6.467 . . . 2 2 2 2
94. Computing matrices dumas.io/PML Numpy’s matrix algebra + recursive splitting to compute ρ(x): import numpy as np def representation(gen_matrices): def _rho(x): if len(x) == 0: # identity matrix return np.eye(2) elif len(x) == 1: # generator matrix return gen_matrices[x] else: N = len(x) return np.dot(_rho(x[:N//2]),_rho(x[N//2:])) return _rho rho = representation( {'a': np.array( [[0,1],[1,1]] )} ) # etc print(rho('aaaaa')) # -> [[3 5], [5 8]]
95. Length and dlength dumas.io/PML eps = 0.0001 rho0 = representation( {'a': np.array( [[0,1],[1,1]] )} ) rho1 = representation( {'a': np.array( [[0,1],[1+eps,1]] )} ) L0 = 2.0*np.arccosh(0.5*np.trace(rho0('aaaaa'))) L1 = 2.0*np.arccosh(0.5*np.trace(rho1('aaaaa'))) dL = (L1 - L0) / eps print(dL) Four calculations like this give the four coordinates of the vector dL associated to a curve.
96. Length and dlength dumas.io/PML eps = 0.0001 rho0 = representation( {'a': np.array( [[0,1],[1,1]] )} ) rho1 = representation( {'a': np.array( [[0,1],[1+eps,1]] )} ) L0 = 2.0*np.arccosh(0.5*np.trace(rho0('aaaaa'))) L1 = 2.0*np.arccosh(0.5*np.trace(rho1('aaaaa'))) dL = (L1 - L0) / eps print(dL) Four calculations like this give the four coordinates of the vector dL associated to a curve. Next problem: 4D space is difficult to visualize!
97. Stereographic projection dumas.io/PML “Open up” the surface in 4-space and flatten to 3-space.
98. Stereographic projection dumas.io/PML “Open up” the surface in 4-space and flatten to 3-space. (p1 , p2 , p3 , p4 ) −→ ( p1 1 − p4 , p2 1 − p4 , p3 1 − p4 )
99. Stereographic projection dumas.io/PML “Open up” the surface in 4-space and flatten to 3-space. (p1 , p2 , p3 , p4 ) −→ ( p1 1 − p4 , p2 1 − p4 , p3 1 − p4 )
100. Stereographic projection dumas.io/PML “Open up” the surface in 4-space and flatten to 3-space. (p1 , p2 , p3 , p4 ) −→ ( p1 1 − p4 , p2 1 − p4 , p3 1 − p4 )
101. Stereographic projection dumas.io/PML “Open up” the surface in 4-space and flatten to 3-space. (p1 , p2 , p3 , p4 ) −→ ( p1 1 − p4 , p2 1 − p4 , p3 1 − p4 )
102. Stereographic projection dumas.io/PML “Open up” the surface in 4-space and flatten to 3-space. (p1 , p2 , p3 , p4 ) −→ ( p1 1 − p4 , p2 1 − p4 , p3 1 − p4 )
103. Rendering dumas.io/PML Projecting each dL to a 3-vector (x, y, z), we obtain data like word L x y z aDBdbd 19.195 0.630 0.729 0.695 aDbcBd 16.790 0.596 0.470 1.133 bcbCBd 14.664 0.177 0.544 1.242 with > 106 rows. To make short curves more prominent, but all curves visible, we draw a sphere centered at (x, y, z) of radius 1/L for each curve.
104. Rendering dumas.io/PML Projecting each dL to a 3-vector (x, y, z), we obtain data like word L x y z aDBdbd 19.195 0.630 0.729 0.695 aDbcBd 16.790 0.596 0.470 1.133 bcbCBd 14.664 0.177 0.544 1.242 with > 106 rows. To make short curves more prominent, but all curves visible, we draw a sphere centered at (x, y, z) of radius 1/L for each curve. Easy way to draw millions of spheres if real-time rendering is not essential?
105. Rendering dumas.io/PML Party like it’s 1992
106. Rendering dumas.io/PML Party like it’s 1992 with ray-tracing!
107. POV-Ray dumas.io/PML POV-Ray is an open-source ray tracer that uses a scene description language (SDL) with syntax reminiscent of C. We generate SDL sphere primitives from the dL data using: outfile.write('sphere { <%f, %f, %f>, %f }\n' % (x,y,z,1.0/L)) → Appending some camera and lighting setup boilerplate gives a full scene file for rendering with POV-Ray.
108. dumas.io/PML Fractal dust pmls05-001
109. dumas.io/PML Rotating the pole (wide angle) pmls05-010
110. dumas.io/PML Rotating the pole (close-up) pmls05-020
111. dumas.io/PML Clifford flow pmls05-030
112. Rings dumas.io/PML
113. Rings dumas.io/PML
114. dumas.io/PML pmls05-071
115. dumas.io/PML Rotating the pole pmls05-041
116. dumas.io/PML Rotating the pole II pmls05-061
117. Twists dumas.io/PML
118. Twists dumas.io/PML
119. Twists dumas.io/PML
120. Twists dumas.io/PML
121. Twists dumas.io/PML
122. Twists dumas.io/PML
123. dumas.io/PML pmls05-081
124. Twists dumas.io/PML
125. Rings poster dumas.io/PML
126. David Dumas david@dumas.io http://dumas.io github.com/daviddumas