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

Lisp as Renaissance Workshop: A Lispy Tour through Mathematical Physics

Lisp as Renaissance Workshop: A Lispy Tour through Mathematical Physics

These slides are from my keynote at the 2022 European Lisp Symposium in Porto, Portugal: https://european-lisp-symposium.org/2022/index.html

The source code for these slides and the demos in the talk live at https://github.com/sritchie/programming-2022.

Abstract:
Lisp is an exquisite medium for the communication of computational ideas. From our most accurate observations of physical reality up through chemistry, biology, and cognition, the universe seems to be computing itself; modeling and simulating these systems in machines has led to incredible technological wealth.

Deep principles and beautiful abstractions seem to drive these systems, but they have always been hard to discover; and we are floundering at the computational frontiers of intelligence, synthetic biology and control systems for our climate. The only way to push forward is to build powerful tools that can communicate and teach.

This talk will take a tour through SICMUtils, a Lisp system designed as a workshop for conducting serious work in mathematical physics and sharing those explorations in a deeply interactive, multiplayer way. The library’s growth parallels our human scientific history; hopefully tools like this will help us write the next chapter.

Cd378611a91eb7852ae19cd582de718a?s=128

Sam Ritchie

March 23, 2022
Tweet

More Decks by Sam Ritchie

Other Decks in Programming

Transcript

  1. LISP AS RENAISSANCE WORKSHOP A Lispy Tour through Mathematical Physics

    , Mentat Collective Sam Ritchie 2
  2. SICMUTILS 4

  3. [(+ (sin 'x) 'x) (+ (sin 12) 2)] [(+ (sin

    x) x) 1.4634270819995652] (tex$$ (up (square (cos (* 't 'phi))) (floor (* 4 'zeta)))) ⎛ ⎝ cos 2 (ϕ t) ⌊4 ζ⌋ ⎞ ⎠ 5
  4. (let [f (literal-function 'f (-> (UP Real Real) Real))] (tex$$

    ((D f) (up 'alpha_1 'alpha_2)))) ⎡ ⎣ ∂0 f ∂1 f ⎛ ⎝ ⎛ ⎝ α 1 α 2 ⎞ ⎠ ⎞ ⎠ ⎛ ⎝ ⎛ ⎝ α 1 α 2 ⎞ ⎠ ⎞ ⎠ ⎤ ⎦ 6
  5. ❤️ OPEN SOURCE ❤️ https://github.com/sicmutils/sicmutils

  6. https://github.com/sritchie/programming-2022 7

  7. AGENDA Why build another CAS? Code as Communication? Literate Systems

    (and their problems) Why it matters What to do? + SICMUtils 9
  8. SCMUTILS BY GJS 11

  9. "programs must be written for people to read, and only

    incidentally for machines to execute." ~ Hal Abelson, Structure and Interpretation of Computer Programs 12
  10. SICM AND FDG 14

  11. MOTIVATION, GOOGLE X 16

  12. DISAPPOINTMENT @ X! 17

  13. 19TH CENTURY SCIENTIFIC COMMUNICATION 19

  14. HISTORY OF VECTOR ANALYSIS 20

  15. 21

  16. CODE AS COMMUNICATION 23

  17. NUMERICAL CODE AS COMMUNICATION? "We personally like Brent's algorithm for

    univariate minimization, as found on pages 79-80 of his book 'Algorithms for Minimization Without Derivatives'. It is pretty reliable and pretty fast, but we cannot explain how it works." ~ scmutils, refman.txt 24
  18. BRENT'S BOOK, 1973 25

  19. FORTRAN: NUMERICAL RECIPES, 1986 26

  20. C++: BOOST, 2006 template <class F, class T> 1 std::pair<T,

    T> brent_find_minima(F f, T min, T max, int bits 2 noexcept(BOOST_MATH_IS_FLOAT(T) && noexcept(std::declval<F 3 { 4 BOOST_MATH_STD_USING 5 bits = (std::min)(policies::digits<T, policies::policy<> > 6 T tolerance = static_cast<T>(ldexp(1.0, 1-bits)); 7 T x; // minima so far 8 T w; // second best point 9 T v; // previous value of w 10 T u; // most recent evaluation point 11 T delta; // The distance moved in the last step 12 T delta2; // The distance moved in the step before last 13 T fu, fv, fw, fx; // function evaluations at u, v, w, x 14 T mid; // midpoint of min and max 15 T fract1 fract2; // minimal relative movement in x 16 27
  21. C++: BOOST, 2006 template <class F, class T> 1 std::pair<T,

    T> brent_find_minima(F f, T min, T max, int bits 2 noexcept(BOOST_MATH_IS_FLOAT(T) && noexcept(std::declval<F 3 { 4 BOOST_MATH_STD_USING 5 bits = (std::min)(policies::digits<T, policies::policy<> > 6 T tolerance = static_cast<T>(ldexp(1.0, 1-bits)); 7 T x; // minima so far 8 T w; // second best point 9 T v; // previous value of w 10 T u; // most recent evaluation point 11 T delta; // The distance moved in the last step 12 T delta2; // The distance moved in the step before last 13 T fu, fv, fw, fx; // function evaluations at u, v, w, x 14 T mid; // midpoint of min and max 15 T fract1 fract2; // minimal relative movement in x 16 } 62 else 63 { 64 // golden section: 65 delta2 = (x >= mid) ? min - x : max - x; 66 delta = golden * delta2; 67 } 68 // update current position: 69 u = (fabs(delta) >= fract1) ? T(x + delta) : (delta > 0 70 fu = f(u); 71 if(fu <= fx) 72 { 73 // good new point is an improvement! 74 // update brackets: 75 if(u >= x) 76 i 77 27
  22. C++: BOOST, 2006 template <class F, class T> 1 std::pair<T,

    T> brent_find_minima(F f, T min, T max, int bits 2 noexcept(BOOST_MATH_IS_FLOAT(T) && noexcept(std::declval<F 3 { 4 BOOST_MATH_STD_USING 5 bits = (std::min)(policies::digits<T, policies::policy<> > 6 T tolerance = static_cast<T>(ldexp(1.0, 1-bits)); 7 T x; // minima so far 8 T w; // second best point 9 T v; // previous value of w 10 T u; // most recent evaluation point 11 T delta; // The distance moved in the last step 12 T delta2; // The distance moved in the step before last 13 T fu, fv, fw, fx; // function evaluations at u, v, w, x 14 T mid; // midpoint of min and max 15 T fract1 fract2; // minimal relative movement in x 16 template <class F, class T> 1 std::pair<T, T> brent_find_minima(F f, T min, T max, int bits 2 noexcept(BOOST_MATH_IS_FLOAT(T) && noexcept(std::declval<F 3 { 4 BOOST_MATH_STD_USING 5 bits = (std::min)(policies::digits<T, policies::policy<> > 6 T tolerance = static_cast<T>(ldexp(1.0, 1-bits)); 7 T x; // minima so far 8 T w; // second best point 9 T v; // previous value of w 10 T u; // most recent evaluation point 11 T delta; // The distance moved in the last step 12 T delta2; // The distance moved in the step before last 13 T fu, fv, fw, fx; // function evaluations at u, v, w, x 14 T mid; // midpoint of min and max 15 T fract1 fract2; // minimal relative movement in x 16 ; x = u; 83 fv = fw; 84 fw = fx; 85 fx = fu; 86 } 87 else 88 { 89 // Oh dear, point u is worse than what we have alrea 90 // even so it *must* be better than one of our endpo 91 if(u < x) 92 min = u; 93 else 94 max = u; 95 if((fu <= fw) || (w == x)) 96 { 97 // however it is at least second best: 98 27
  23. C++: BOOST, 2006 template <class F, class T> 1 std::pair<T,

    T> brent_find_minima(F f, T min, T max, int bits 2 noexcept(BOOST_MATH_IS_FLOAT(T) && noexcept(std::declval<F 3 { 4 BOOST_MATH_STD_USING 5 bits = (std::min)(policies::digits<T, policies::policy<> > 6 T tolerance = static_cast<T>(ldexp(1.0, 1-bits)); 7 T x; // minima so far 8 T w; // second best point 9 T v; // previous value of w 10 T u; // most recent evaluation point 11 T delta; // The distance moved in the last step 12 T delta2; // The distance moved in the step before last 13 T fu, fv, fw, fx; // function evaluations at u, v, w, x 14 T mid; // midpoint of min and max 15 T fract1 fract2; // minimal relative movement in x 16 template <class F, class T> 1 std::pair<T, T> brent_find_minima(F f, T min, T max, int bits 2 noexcept(BOOST_MATH_IS_FLOAT(T) && noexcept(std::declval<F 3 { 4 BOOST_MATH_STD_USING 5 bits = (std::min)(policies::digits<T, policies::policy<> > 6 T tolerance = static_cast<T>(ldexp(1.0, 1-bits)); 7 T x; // minima so far 8 T w; // second best point 9 T v; // previous value of w 10 T u; // most recent evaluation point 11 T delta; // The distance moved in the last step 12 T delta2; // The distance moved in the step before last 13 T fu, fv, fw, fx; // function evaluations at u, v, w, x 14 T mid; // midpoint of min and max 15 T fract1 fract2; // minimal relative movement in x 16 template <class F, class T> 1 std::pair<T, T> brent_find_minima(F f, T min, T max, int bits 2 noexcept(BOOST_MATH_IS_FLOAT(T) && noexcept(std::declval<F 3 { 4 BOOST_MATH_STD_USING 5 bits = (std::min)(policies::digits<T, policies::policy<> > 6 T tolerance = static_cast<T>(ldexp(1.0, 1-bits)); 7 T x; // minima so far 8 T w; // second best point 9 T v; // previous value of w 10 T u; // most recent evaluation point 11 T delta; // The distance moved in the last step 12 T delta2; // The distance moved in the step before last 13 T fu, fv, fw, fx; // function evaluations at u, v, w, x 14 T mid; // midpoint of min and max 15 T fract1 fract2; // minimal relative movement in x 16 w = u; 100 fv = fw; 101 fw = fu; 102 } 103 else if((fu <= fv) || (v == x) || (v == w)) 104 { 105 // third best: 106 v = u; 107 fv = fu; 108 } 109 } 110 111 }while(--count); 112 113 max iter -= count; 114 27
  24. C++: BOOST, 2006 template <class F, class T> 1 std::pair<T,

    T> brent_find_minima(F f, T min, T max, int bits 2 noexcept(BOOST_MATH_IS_FLOAT(T) && noexcept(std::declval<F 3 { 4 BOOST_MATH_STD_USING 5 bits = (std::min)(policies::digits<T, policies::policy<> > 6 T tolerance = static_cast<T>(ldexp(1.0, 1-bits)); 7 T x; // minima so far 8 T w; // second best point 9 T v; // previous value of w 10 T u; // most recent evaluation point 11 T delta; // The distance moved in the last step 12 T delta2; // The distance moved in the step before last 13 T fu, fv, fw, fx; // function evaluations at u, v, w, x 14 T mid; // midpoint of min and max 15 T fract1 fract2; // minimal relative movement in x 16 template <class F, class T> 1 std::pair<T, T> brent_find_minima(F f, T min, T max, int bits 2 noexcept(BOOST_MATH_IS_FLOAT(T) && noexcept(std::declval<F 3 { 4 BOOST_MATH_STD_USING 5 bits = (std::min)(policies::digits<T, policies::policy<> > 6 T tolerance = static_cast<T>(ldexp(1.0, 1-bits)); 7 T x; // minima so far 8 T w; // second best point 9 T v; // previous value of w 10 T u; // most recent evaluation point 11 T delta; // The distance moved in the last step 12 T delta2; // The distance moved in the step before last 13 T fu, fv, fw, fx; // function evaluations at u, v, w, x 14 T mid; // midpoint of min and max 15 T fract1 fract2; // minimal relative movement in x 16 template <class F, class T> 1 std::pair<T, T> brent_find_minima(F f, T min, T max, int bits 2 noexcept(BOOST_MATH_IS_FLOAT(T) && noexcept(std::declval<F 3 { 4 BOOST_MATH_STD_USING 5 bits = (std::min)(policies::digits<T, policies::policy<> > 6 T tolerance = static_cast<T>(ldexp(1.0, 1-bits)); 7 T x; // minima so far 8 T w; // second best point 9 T v; // previous value of w 10 T u; // most recent evaluation point 11 T delta; // The distance moved in the last step 12 T delta2; // The distance moved in the step before last 13 T fu, fv, fw, fx; // function evaluations at u, v, w, x 14 T mid; // midpoint of min and max 15 T fract1 fract2; // minimal relative movement in x 16 template <class F, class T> 1 std::pair<T, T> brent_find_minima(F f, T min, T max, int bits 2 noexcept(BOOST_MATH_IS_FLOAT(T) && noexcept(std::declval<F 3 { 4 BOOST_MATH_STD_USING 5 bits = (std::min)(policies::digits<T, policies::policy<> > 6 T tolerance = static_cast<T>(ldexp(1.0, 1-bits)); 7 T x; // minima so far 8 T w; // second best point 9 T v; // previous value of w 10 T u; // most recent evaluation point 11 T delta; // The distance moved in the last step 12 T delta2; // The distance moved in the step before last 13 T fu, fv, fw, fx; // function evaluations at u, v, w, x 14 T mid; // midpoint of min and max 15 T fract1 fract2; // minimal relative movement in x 16 fw = fu; 102 } 103 else if((fu <= fv) || (v == x) || (v == w)) 104 { 105 // third best: 106 v = u; 107 fv = fu; 108 } 109 } 110 111 }while(--count); 112 113 max_iter -= count; 114 115 return std::make_pair(x, fx); 116 } 117 27
  25. PYTHON: SCIPY, 2001 def optimize(self): 1 # set up for

    optimization 2 func = self.func 3 xa, xb, xc, fa, fb, fc, funcalls = self.get_bracket_info() 4 _mintol = self._mintol 5 _cg = self._cg 6 ################################# 7 #BEGIN CORE ALGORITHM 8 ################################# 9 x = w = v = xb 10 fw = fv = fx = func(*((x,) + self.args)) 11 if (xa < xc): 12 a = xa 13 b = xc 14 else: 15 a = xc 16 28
  26. SCHEME: 1987 ;;; Brent's algorithm for univariate minimization -- transcrib

    1 ;;; pages 79-80 of his book "Algorithms for Minimization Witho 2 3 (define (brent-min f a b eps) 4 (let ((a (min a b)) (b (max a b)) 5 (maxcount 100) 6 (small-bugger-factor *sqrt-machine-epsilon*) 7 (g (/ (- 3 (sqrt 5)) 2)) 8 (d 0) (e 0) (old-e 0) (p 0) (q 0) (u 0) (fu 0)) 9 (let* ((x (+ a (* g (- b a)))) 10 (fx (f x)) 11 (w x) (fw fx) (v x) (fv fx)) 12 (let loop ((count 0)) 13 (if (> count maxcount) 14 (list 'maxcount x fx count) ;failed to converge 15 (let* ((tol (+ (* eps (abs x)) small-bugger-factor)) 16 29
  27. ACTUAL CORE IDEA 30

  28. SO WHAT? What's wrong with this? How can we do

    better? 31
  29. EXISTING IDEAS Literate Programming LOGO's Microworlds Notebooks, Mathematica 33

  30. 34

  31. LOGO'S MICROWORLDS 35

  32. 36

  33. NOTEBOOKS, MATHEMATICA 37

  34. WHY AREN'T THESE WORKING? Literate Programming == one-way Science is

    Multiplayer Dynamic is Too Seductive Real Work doesn't happen here 39
  35. WHY AREN'T THESE WORKING? Literate Programming == one-way Science is

    Multiplayer Dynamic is Too Seductive Real Work doesn't happen here 39
  36. WHY AREN'T THESE WORKING? Literate Programming == one-way Science is

    Multiplayer Dynamic is Too Seductive Real Work doesn't happen here 39
  37. WHY AREN'T THESE WORKING? Literate Programming == one-way Science is

    Multiplayer Dynamic is Too Seductive Real Work doesn't happen here 39
  38. WHY AREN'T THESE WORKING? Literate Programming == one-way Science is

    Multiplayer Dynamic is Too Seductive Real Work doesn't happen here 39
  39. WHY CARE? 41

  40. "FUTURE OF EDUCATION" Mythic Understanding Romantic Understanding Philosophic Understanding Ironic

    Understanding 42
  41. WHAT TO DO? 44

  42. SICM AND FDG AS CLUES 45

  43. 46

  44. EXPLANATION AS SIDE EFFECT 47

  45. EULER-LAGRANGE EQUATIONS S[q](t a , t b ) = D

    ∫ t b t a L(t, q(t), Dq(t))dt D ∫ t b t a L(t, q(t), Dq(t))dt = 0 49
  46. EULER-LAGRANGE EQUATIONS d dt ∂L ∂ ˙ q − ∂L

    ∂q = 0 50
  47. EULER-LAGRANGE EQUATIONS "What could this expression possibly mean?" d dt

    ∂L ∂ ˙ q − ∂L ∂q = 0 50
  48. EXPAND: d dt ∂L ∂ ˙ q − ∂L ∂q

    = 0 51
  49. EXPAND: d dt ∂L ∂ ˙ q − ∂L ∂q

    = 0 d dt ∂L(t, q, ˙ q) ∂ ˙ q − ∂L(t, q, ˙ q) ∂q = 0 ⎛ ⎝ q=w(t) ˙ q= dw(t) dt ⎞ ⎠ q=w(t) ˙ q= dw(t) dt 51
  50. OKAY, FINE d dt ∂L ∂ ˙ q − ∂L

    ∂q = 0 d dt ((∂ 2 L) (t, w(t), d dt w(t))) −(∂ 1 L) (t, w(t), d dt w(t)) = 0 52
  51. SUBSTITUTIONS (Df)(t) = d dx f(x) x=t Γ[w](t) = (t,

    w(t), d dt w(t)) 53
  52. d dt ((∂ 2 L) (Γ[w](t))) − (∂ 1 L)

    (Γ[w](t)) = 0 54
  53. d dt ((∂ 2 L) (Γ[w](t))) − (∂ 1 L)

    (Γ[w](t)) = 0 D ((∂ 2 L) ∘ (Γ[w])) − (∂ 1 L) ∘ (Γ[w]) = 0 54
  54. d dt ((∂ 2 L) (Γ[w](t))) − (∂ 1 L)

    (Γ[w](t)) = 0 D ((∂ 2 L) ∘ (Γ[w])) − (∂ 1 L) ∘ (Γ[w]) = 0 (defn Lagrange-equations [L] 1 (fn [w] 2 (- (D (comp ((partial 2) L) (Gamma w))) 3 (comp ((partial 1) L) (Gamma w))))) 4 54
  55. d dt ((∂ 2 L) (Γ[w](t))) − (∂ 1 L)

    (Γ[w](t)) = 0 D ((∂ 2 L) ∘ (Γ[w])) − (∂ 1 L) ∘ (Γ[w]) = 0 (defn Lagrange-equations [L] 1 (fn [w] 2 (- (D (comp ((partial 2) L) (Gamma w))) 3 (comp ((partial 1) L) (Gamma w))))) 4 (fn [w] (- (D (comp ((partial 2) L) (Gamma w))) (comp ((partial 1) L) (Gamma w))))) (defn Lagrange-equations [L] 1 2 3 4 54
  56. (let [L (- (comp (literal-function 'T) velocity) (comp (literal-function 'V)

    coordinate)) w (literal-function 'w)] (tex$$ (((Lagrange-equations L) w) 't))) D 2 w (t) D 2 T (Dw (t)) + DV (w (t)) 55
  57. SIMPLE HARMONIC MOTION (- (* 1/2 m (square v)) (defn

    L-harmonic 1 "Returns a Lagrangian of a simple harmonic oscillator (mass-s 2 3 m is the mass and k is the spring constant used in Hooke's la 4 [m k] 5 (fn [[_ q v]] 6 7 (* 1/2 k (square q))))) 8 (defn proposed-solution [t] (* 'a (cos (+ (* 'omega t) 'phi)))) (let [L (L-harmonic 'm 'k) w proposed-solution] (tex$$ (((Lagrange-equations L) w) 't))) −a m ω 2 cos (ω t + ϕ) + a k cos (ω t + ϕ) 56
  58. SIMPLE HARMONIC MOTION (- (* 1/2 m (square v)) (defn

    L-harmonic 1 "Returns a Lagrangian of a simple harmonic oscillator (mass-s 2 3 m is the mass and k is the spring constant used in Hooke's la 4 [m k] 5 (fn [[_ q v]] 6 7 (* 1/2 k (square q))))) 8 (* 1/2 k (square q))))) (defn L-harmonic 1 "Returns a Lagrangian of a simple harmonic oscillator (mass-s 2 3 m is the mass and k is the spring constant used in Hooke's la 4 [m k] 5 (fn [[_ q v]] 6 (- (* 1/2 m (square v)) 7 8 (defn proposed-solution [t] (* 'a (cos (+ (* 'omega t) 'phi)))) (let [L (L-harmonic 'm 'k) w proposed-solution] (tex$$ (((Lagrange-equations L) w) 't))) −a m ω 2 cos (ω t + ϕ) + a k cos (ω t + ϕ) 56
  59. HOOKE'S LAW (let [L (L-harmonic 'm 'k) w (literal-function 'w)]

    (tex$$ (((Lagrange-equations L) w) 't))) k w (t) + m D 2 w (t) 57
  60. CLERK DEMO Let's convince ourselves that that is true, by

    doing the spring demo in the browser. 59
  61. COMMUNITY Sussman et. al. Colin Smith Clerk: Martin Kavalar, Jack

    Rusher, Nextjournal team SCI: Michiel Borkent (@borkdude) Mathbox: Chris Chudzicki, Steven Wittens 61
  62. HOW TO GET INVOLVED? WHAT'S NEXT? SICM and FDG: Executable

    Textbooks Full library as essays Collaborative editing, simulation Please Steal! 63
  63. THANKS! Sam Ritchie, Mentat Collective Slides, Demos live at @sritchie

    https://github.com/sritchie/programming-2022 65
  64. None