Comparison to Scheme scmutils — org.mentat/emmy 0.32.0

Comparison to scmutils

Rather than just quasi-mechanically translate the Scheme to Clojure, We’ve studied the implementation of the system before bringing it to Clojure, and have used TDD throughout the project (which turned out to be absolutely essential as I considered various approaches to problems posed by the Scheme code base). At this writing there are over 4200 unit tests.

Here’s a side-by-side example of scmutils Scheme code and SICMUtils Clojure code. First, the Scheme:

;; Scheme
(define ((L-central-polar m U) local)
  (let ((q (coordinate local))
        (qdot (velocity local)))
    (let ((r (ref q 0)) (phi (ref q 1))
          (rdot (ref qdot 0)) (phidot (ref qdot 1)))
      (- (* 1/2 m
           (+ (square rdot)
              (square (* r phidot))) )
         (U r)))))

Then the same function in idiomatic Clojure:

;; Clojure
(defn L-central-polar [m U]
  (fn [[_ [r] [rdot φdot]]]
    (- (* 1/2 m
          (+ (square rdot)
             (square (* r φdot))))
       (U r))))

We can see a few things from this example. L-central-polar wants to compute a Lagrangian for a point mass m in a potential field U. In Scheme, it’s possible to specify currying at the site of a function’s definition: (L-central-polar m U) returns a function of the local tuple (a sequence of time, generalized coordinates, and generalized velocities). We don’t have that syntax in Clojure, but instead have something even more useful: argument destructuring. We can pick out exactly the coordinates we want out of the local tuple components directly.

While function definitions cannot be typed directly from the book, function applications in Clojure and Scheme are the same. The following works in both systems:

(((Lagrange-equations (L-central-polar 'm (literal-function 'U)))
  (up (literal-function 'r)
      (literal-function 'φ)))
  't)

yielding:

(down
 (+ (* -1N m (expt ((D φ) t) 2) (r t)) (* m (((expt D 2) r) t)) ((D U) (r t)))
 (+ (* 2N m ((D φ) t) (r t) ((D r) t)) (* m (expt (r t) 2) (((expt D 2) φ) t))))

Which, modulo a few things, is what Scmutils would give. From later in SICM (pp. 81-2) we have, in Scheme:

(define ((T3-spherical m) state)
  (let ((t (time state))
        (q (coordinate state))
        (qdot (velocity state)))
    (let ((r (ref q 0))
          (theta (ref q 1))
          (phi (ref q 2))
          (rdot (ref qdot 0))
          (thetadot (ref qdot 1))
          (phidot (ref qdot 2)))
      (* 1/2 m
         (+ (square rdot)
            (square (* r thetadot))
            (square (* r (sin theta) phidot)))))))

(define (L3-central m Vr)
  (define (Vs state)
    (let ((r (ref (coordinate state) 0)))
      (Vr r)))
  (- (T3-spherical m) Vs))

(((partial 1) (L3-central ’m (literal-function ’V)))
  (up ’t
      (up ’r ’theta ’phi)
      (up ’rdot ’thetadot ’phidot)))

And in Clojure, using a couple of simplifying definitions:

(def V (literal-function 'V))
(def spherical-state (up 't
                         (up 'r 'θ 'φ)
                         (up 'rdot 'θdot 'φdot)))
(defn T3-spherical [m]
  (fn [[t [r θ φ] [rdot θdot φdot]]]
    (* 1/2 m (+ (square rdot)
                (square (* r θdot))
                (square (* r (sin θ) φdot))))))

(defn L3-central [m Vr]
  (let [Vs (fn [[_ [r]]] (Vr r))]
    (- (T3-spherical m) Vs)))

(((partial 1) (L3-central 'm V)) spherical-state)

yielding

(down
 (+ (* m r (expt φdot 2) (expt (sin θ) 2)) (* m r (expt θdot 2)) (* -1 ((D V) r)))
 (* m (expt r 2) (expt φdot 2) (sin θ) (cos θ))
 0)

Which again agrees with Scmutils modulo notation. (These results are examples of "down tuples", or covariant vectors, since they represent derivatives of objects in primal space.)

The partial derivative operation is called partial in Scmutils, but Clojure defines partial to mean partial function application. In this system, we take a page from JavaScript and replace partial with a shim which will compute partial derivatives when all the arguments are integers and fall back to Clojure’s definition of partial otherwise. Since it doesn’t make sense to partially apply an integer, partial should just do the right thing.

You could render that result in TeX:

TeX result

using the →TeX function. You can also use →infix to obtain:

down(m r φdot² sin²(θ) + m r θdot² - DV(r), m r² φdot² sin(θ) cos(θ), 0)

or even →JavaScript to get:

function(D, V, m, r, θ, θdot, φdot) {
  var _0001 = Math.sin(θ);
  var _0002 = Math.pow(φdot, 2);
  return [m * r * _0002 * Math.pow(_0001, 2) + m * r * Math.pow(θdot, 2) - (D(V)(r)), m * Math.pow(r, 2) * _0002 * _0001 * Math.cos(θ), 0];
}

(For rendering into code, a simple common-subexpression extraction algorithm is used.)

Numerical Methods

This system uses the delightful Apache Commons Math library to implement some of the numerical methods needed for the book. Many of the methods we used to pull from Commons are now implemented in native Clojure(script); our goal is to remove this dependency once we get a pure Clojure implementation of the Gragg-Bulirsch-Stoer ODE integrator implemented.

The Scmutils simplifier has three engines: a polynomial-based simplifier useful for grouping like terms, a rational-function-based simplifier used for cancellation in fractional expressions, and a rule-based simplifier to apply identities like sin² x + cos² x = 1.

I have implemented all of these, but acceptable performance from the rational-function simplifier is waiting for an implementation of Zippel’s algorithm for fast multivariate polynomial GCD operations. Currently we use a recursive Euclid algorithm, which gives acceptable results for expressions of medium complexity, but there is more to be done.