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
   

   View Slide

2. SICMUTILS
   4
   

   View Slide

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
   

   View Slide

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
   

   View Slide

5. ❤️
   OPEN SOURCE ❤️
   https://github.com/sicmutils/sicmutils
   

   View Slide

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

   View Slide

7. AGENDA
   Why build another CAS?
   Code as Communication?
   Literate Systems (and their problems)
   Why it matters
   What to do? + SICMUtils
   9
   

   View Slide

8. SCMUTILS BY GJS
   11
   

   View Slide

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
   

   View Slide

10. SICM AND FDG
    
    
    14
    

    View Slide

11. MOTIVATION, GOOGLE X
    16
    

    View Slide

12. DISAPPOINTMENT @ X!
    17
    

    View Slide

13. 19TH CENTURY SCIENTIFIC
    COMMUNICATION
    19
    

    View Slide

14. HISTORY OF VECTOR ANALYSIS
    20
    

    View Slide

15. View Slide

16. CODE AS COMMUNICATION
    23
    

    View Slide

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
    

    View Slide

18. BRENT'S BOOK, 1973
    25
    

    View Slide

19. FORTRAN: NUMERICAL RECIPES, 1986
    26
    

    View Slide

20. C++: BOOST, 2006
    template
    1
    std::pair brent_find_minima(F f, T min, T max, int bits
    2
    noexcept(BOOST_MATH_IS_FLOAT(T) && noexcept(std::declval3
    {
    4
    BOOST_MATH_STD_USING
    5
    bits = (std::min)(policies::digits >
    6
    T tolerance = static_cast(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
    

    View Slide

21. C++: BOOST, 2006
    template
    1
    std::pair brent_find_minima(F f, T min, T max, int bits
    2
    noexcept(BOOST_MATH_IS_FLOAT(T) && noexcept(std::declval3
    {
    4
    BOOST_MATH_STD_USING
    5
    bits = (std::min)(policies::digits >
    6
    T tolerance = static_cast(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
    

    View Slide

22. C++: BOOST, 2006
    template
    1
    std::pair brent_find_minima(F f, T min, T max, int bits
    2
    noexcept(BOOST_MATH_IS_FLOAT(T) && noexcept(std::declval3
    {
    4
    BOOST_MATH_STD_USING
    5
    bits = (std::min)(policies::digits >
    6
    T tolerance = static_cast(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
    1
    std::pair brent_find_minima(F f, T min, T max, int bits
    2
    noexcept(BOOST_MATH_IS_FLOAT(T) && noexcept(std::declval3
    {
    4
    BOOST_MATH_STD_USING
    5
    bits = (std::min)(policies::digits >
    6
    T tolerance = static_cast(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
    

    View Slide

23. C++: BOOST, 2006
    template
    1
    std::pair brent_find_minima(F f, T min, T max, int bits
    2
    noexcept(BOOST_MATH_IS_FLOAT(T) && noexcept(std::declval3
    {
    4
    BOOST_MATH_STD_USING
    5
    bits = (std::min)(policies::digits >
    6
    T tolerance = static_cast(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
    1
    std::pair brent_find_minima(F f, T min, T max, int bits
    2
    noexcept(BOOST_MATH_IS_FLOAT(T) && noexcept(std::declval3
    {
    4
    BOOST_MATH_STD_USING
    5
    bits = (std::min)(policies::digits >
    6
    T tolerance = static_cast(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
    1
    std::pair brent_find_minima(F f, T min, T max, int bits
    2
    noexcept(BOOST_MATH_IS_FLOAT(T) && noexcept(std::declval3
    {
    4
    BOOST_MATH_STD_USING
    5
    bits = (std::min)(policies::digits >
    6
    T tolerance = static_cast(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
    

    View Slide

24. C++: BOOST, 2006
    template
    1
    std::pair brent_find_minima(F f, T min, T max, int bits
    2
    noexcept(BOOST_MATH_IS_FLOAT(T) && noexcept(std::declval3
    {
    4
    BOOST_MATH_STD_USING
    5
    bits = (std::min)(policies::digits >
    6
    T tolerance = static_cast(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
    1
    std::pair brent_find_minima(F f, T min, T max, int bits
    2
    noexcept(BOOST_MATH_IS_FLOAT(T) && noexcept(std::declval3
    {
    4
    BOOST_MATH_STD_USING
    5
    bits = (std::min)(policies::digits >
    6
    T tolerance = static_cast(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
    1
    std::pair brent_find_minima(F f, T min, T max, int bits
    2
    noexcept(BOOST_MATH_IS_FLOAT(T) && noexcept(std::declval3
    {
    4
    BOOST_MATH_STD_USING
    5
    bits = (std::min)(policies::digits >
    6
    T tolerance = static_cast(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
    1
    std::pair brent_find_minima(F f, T min, T max, int bits
    2
    noexcept(BOOST_MATH_IS_FLOAT(T) && noexcept(std::declval3
    {
    4
    BOOST_MATH_STD_USING
    5
    bits = (std::min)(policies::digits >
    6
    T tolerance = static_cast(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
    

    View Slide

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
    

    View Slide

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
    

    View Slide

27. ACTUAL CORE IDEA
    30
    

    View Slide

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

    View Slide

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

    View Slide

30. View Slide

31. LOGO'S MICROWORLDS
    35
    

    View Slide

32. View Slide

33. NOTEBOOKS, MATHEMATICA
    
    
    37
    

    View Slide

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

    View Slide

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

    View Slide

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

    View Slide

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

    View Slide

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

    View Slide

39. WHY CARE?
    41
    

    View Slide

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

    View Slide

41. WHAT TO DO?
    44
    

    View Slide

42. SICM AND FDG AS CLUES
    
    
    45
    

    View Slide

43. View Slide

44. EXPLANATION AS SIDE EFFECT
    47
    

    View Slide

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
    

    View Slide

46. EULER-LAGRANGE EQUATIONS
    d
    dt
    ∂L
    ∂ ˙
    q
    −
    ∂L
    ∂q
    = 0
    50
    

    View Slide

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

    View Slide

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

    View Slide

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
    

    View Slide

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
    

    View Slide

51. SUBSTITUTIONS
    (Df)(t) =
    d
    dx
    f(x)
    x=t
    Γ[w](t) = (t, w(t),
    d
    dt
    w(t))
    53
    

    View Slide

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

    View Slide

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

    View Slide

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
    

    View Slide

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
    

    View Slide

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
    

    View Slide

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
    

    View Slide

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
    

    View Slide

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
    

    View Slide

60. CLERK DEMO
    Let's convince ourselves that that is true, by doing the
    spring
    demo in the browser.
    59
    

    View Slide

61. COMMUNITY
    Sussman et. al.
    Colin Smith
    Clerk: Martin Kavalar, Jack Rusher, Nextjournal
    team
    SCI: Michiel Borkent (@borkdude)
    Mathbox: Chris Chudzicki, Steven Wittens
    61
    

    View Slide

62. HOW TO GET INVOLVED? WHAT'S NEXT?
    SICM and FDG: Executable Textbooks
    Full library as essays
    Collaborative editing, simulation
    Please Steal!
    63
    

    View Slide

63. THANKS!
    Sam Ritchie, Mentat Collective
    Slides, Demos live at
    @sritchie
    https://github.com/sritchie/programming-2022
    65
    

    View Slide

64. View Slide