Upgrade to Pro — share decks privately, control downloads, hide ads and more …

David Dumas - Python for mathematical visualiza...

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/

PyCon 2017

May 21, 2017
Tweet

More Decks by PyCon 2017

Other Decks in Programming

Transcript

  1. Goal Describe a mathematical visualization project that uses Python. Highlight

    characteristics of Python that make it suitable for this application.
  2. 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))
  3. 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))
  4. 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()
  5. 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
  6. 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
  7. Twisting is search & replace dumas.io/PML Tcd is equivalent to

    substituting: a → a b → b c → d d → Dcd e → e
  8. 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
  9. 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.
  10. 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)
  11. 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.
  12. 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.
  13. 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]]
  14. 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.
  15. 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!
  16. 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 )
  17. 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 )
  18. 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 )
  19. 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 )
  20. 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 )
  21. 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.
  22. 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?
  23. 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.