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

Lisp as Renaissance Workshop: A Lispy Tour thro...

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.

Sam Ritchie

March 23, 2022
Tweet

More Decks by Sam Ritchie

Other Decks in Programming

Transcript

  1. [(+ (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
  2. (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
  3. AGENDA Why build another CAS? Code as Communication? Literate Systems

    (and their problems) Why it matters What to do? + SICMUtils 9
  4. "programs must be written for people to read, and only

    incidentally for machines to execute." ~ Hal Abelson, Structure and Interpretation of Computer Programs 12
  5. 21

  6. 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
  7. 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
  8. 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
  9. 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
  10. 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
  11. 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
  12. 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
  13. 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
  14. 34

  15. 36

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

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

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

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

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

    Multiplayer Dynamic is Too Seductive Real Work doesn't happen here 39
  21. 46

  22. 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
  23. 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
  24. 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
  25. d dt ((∂ 2 L) (Γ[w](t))) − (∂ 1 L)

    (Γ[w](t)) = 0 D ((∂ 2 L) ∘ (Γ[w])) − (∂ 1 L) ∘ (Γ[w]) = 0 54
  26. 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
  27. 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
  28. (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
  29. 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
  30. 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
  31. 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
  32. CLERK DEMO Let's convince ourselves that that is true, by

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

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

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

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