110

# 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.

March 23, 2022

## Transcript

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

, Mentat Collective Sam Ritchie 2

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

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

(and their problems) Why it matters What to do? + SICMUtils 9

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

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

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

better? 31

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

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

Understanding 42

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

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

∂L ∂ ˙ q − ∂L ∂q = 0 50

= 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