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 full-size slide

  2. [(+ (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 full-size slide

  3. (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 full-size slide

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

    View full-size slide

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

    View full-size slide

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

    View full-size slide

  7. SCMUTILS BY GJS
    11

    View full-size slide

  8. "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 full-size slide

  9. SICM AND FDG


    14

    View full-size slide

  10. MOTIVATION, GOOGLE X
    16

    View full-size slide

  11. DISAPPOINTMENT @ X!
    17

    View full-size slide

  12. 19TH CENTURY SCIENTIFIC
    COMMUNICATION
    19

    View full-size slide

  13. HISTORY OF VECTOR ANALYSIS
    20

    View full-size slide

  14. CODE AS COMMUNICATION
    23

    View full-size slide

  15. 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 full-size slide

  16. BRENT'S BOOK, 1973
    25

    View full-size slide

  17. FORTRAN: NUMERICAL RECIPES, 1986
    26

    View full-size slide

  18. 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 full-size slide

  19. 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 full-size 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
    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 full-size 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
    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 full-size 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
    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 full-size slide

  23. 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 full-size slide

  24. 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 full-size slide

  25. ACTUAL CORE IDEA
    30

    View full-size slide

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

    View full-size slide

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

    View full-size slide

  28. LOGO'S MICROWORLDS
    35

    View full-size slide

  29. NOTEBOOKS, MATHEMATICA


    37

    View full-size slide

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

    View full-size slide

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

    View full-size slide

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

    View full-size slide

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

    View full-size 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 full-size slide

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

    View full-size slide

  36. WHAT TO DO?
    44

    View full-size slide

  37. SICM AND FDG AS CLUES


    45

    View full-size slide

  38. EXPLANATION AS SIDE EFFECT
    47

    View full-size slide

  39. 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 full-size slide

  40. EULER-LAGRANGE EQUATIONS
    d
    dt
    ∂L
    ∂ ˙
    q

    ∂L
    ∂q
    = 0
    50

    View full-size slide

  41. EULER-LAGRANGE EQUATIONS
    "What could this expression possibly mean?"
    d
    dt
    ∂L
    ∂ ˙
    q

    ∂L
    ∂q
    = 0
    50

    View full-size slide

  42. EXPAND:
    d
    dt
    ∂L
    ∂ ˙
    q

    ∂L
    ∂q
    = 0
    51

    View full-size slide

  43. 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 full-size slide

  44. 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 full-size slide

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

    View full-size slide

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

    View full-size slide

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

    View full-size slide

  48. 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 full-size slide

  49. 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 full-size slide

  50. (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 full-size slide

  51. 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 full-size slide

  52. 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 full-size slide

  53. 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 full-size slide

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

    View full-size slide

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

    View full-size slide

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

    View full-size slide

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

    View full-size slide