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.

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

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

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

    View Slide

  31. LOGO'S MICROWORLDS
    35

    View Slide

  32. 36

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

    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