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.

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
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
31. LOGO'S MICROWORLDS 35
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
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