Implement unit systems and SI Units, fold in constants from scmutils

[teodorlu]

Implement unit systems and SI Units, fold in constants from scmutils

[sritchie]

Okay, I made a bunch of notes on how to do this in Slack before realizing I should copy them over. Here we go!

@teodorlu, for when you’re around, I’m out of code mode for real this time!
you’ve made me realize that a tour of how to navigate the scheme library is in order for anyone interested.
a sketch here - I dockerized the mechanics library, ie, scmutils, and it lives here: https://hub.docker.com/r/sritchie/mechanics
If you have docker installed, make this file executable and put it somewhere on your path:

#!/bin/bash
if xhost >& /dev/null ; then
  # NOTE - On OS X, you have to have enabled network connections in XQuartz!
  xhost + 127.0.0.1
fi
workdir="$PWD"
docker run \
       --ipc host \
       --interactive --tty --rm \
       --workdir $workdir \
       --volume $workdir:$workdir \
       -e DISPLAY=host.docker.internal:0 \
       sritchie/mechanics "$@"

the docker page I linked has some instructions if you want to get latex expressions showing or animations working, but I VERY often am opening up this program and fiddling, testing the scmutils implementation by using it directly

next, I have a local copy of the source (here is the 2015 version, ~almost totally the same https://github.com/Tipoca/scmutils/tree/master/src)

“kernel” is where most of the “basic”, bread and butter code lives - units lives in its own folder, but I have a rough mental map of where he calls it (edited)

there is no good module system so I am constantly doing “find in project” in emacs, with projectile. If you search for “define” before a function naem that jumps to impl fairly well (I’m sure I could get some CTAGS thing going but I’ve been fumbling along without it)

then we start digging into a promising file like with-units.scm… https://github.com/Tipoca/scmutils/blob/master/src/units/with-units.scm

these are the money statements:

(assign-operation 'sin              u:sin              angular?)
(assign-operation 'cos              u:cos              angular?)
(assign-operation 'exp              u:cos              angular?)

assign-operation is where the generic operations get installed. So here we see that there is some predicate, angular?, and these u:sin etc functions know how to process them

lines like this:

(assign-operation '*   u:*u    with-units?               units?)
(assign-operation '*   u:u*    units?                    with-units?)

show that this functions, u:*u etc multiply a thing with units by a “unit” object

then he sticks these little examples all over as “tests”

(pe (definite-integral
      (lambda (r)
	(/ (* :G earth-mass (& 1 &kilogram))
	   (square (+ earth-radius r))))
      (& 0 &meter) (& 1 &meter)))
(& 9.824031599863007 &joule)

that & thing is a function:

;;; & is used to attach units to a number, or to check that a number
;;; has the given units.
(define (& value u1 #!optional u2)
  (let ((units (if (default-object? u2) u1 u2))
        (scale (if (default-object? u2) 1 u1)))
    (assert (and (not (units? value)) (number? scale) (units? units)))
    (if (with-units? value)
        (if (equal? (unit-exponents units)
                    (unit-exponents (u:units value)))
            value
            (error "Units do not match: &" value units))
        (with-units (g:* scale (unit-scale units) value)
                    (make-unit (unit-system units)
                               (unit-exponents units)
                               1)))))
(define *unit-constructor* '&)

and then you kind of squint, reformat, pretend its clojure and start going!

(with-units (g:* scale (unit-scale units) value)
            (make-unit (unit-system units)
                       (unit-exponents units)
                       1))

@teodorlu so looks like the functions all referenced here are key. we have a system, exponents, scale, with-units and make-unit.

[littleredcomputer]

Sam's recent introduction of a metadata slot for expressions might be useful here: could we include a slot in the metadata for units? We can use the simplifier we have to canonicalize units to their L T M (etc) form. (in fact our simplifier is overkill because we don't have addition to worry about.

…

On Mon, Dec 7, 2020 at 1:47 PM Sam Ritchie ***@***.***> wrote: @teodorlu <https://github.com/teodorlu> , I'll write you tonight, sorry for the delay, lots going on. I am going to start accumulating some prior art here that we can look at for design, if we want to move beyond / take a different tack than scmutils has taken. JS - https://github.com/sharkdp/insect Clojure - https://github.com/arrdem/meajure - https://github.com/mfey/units2 - https://github.com/fogus/minderbinder - https://github.com/martintrojer/frinj/ — You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub <https://github.com/littleredcomputer/sicmutils/issues/181#issuecomment-740199513>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AAB3Y7MLDHUMFUQDANNF4GTSTVEO5ANCNFSM4UO5I4NA> .

[sritchie]

@teodorlu , I'm assigning this to you, which gives you full license to ask me any questions you like! This is really exciting.

[sritchie]

More prior art, with types, in haskell! http://hackage.haskell.org/package/dimensional

[sritchie]

@teodorlu , good questions! The main thing to note (which we'll look at together) is that & is a function in scmutils:

;;; & is used to attach units to a number, or to check that a number
;;; has the given units.

(define (& value u1 #!optional u2)
  (let ((units (if (default-object? u2) u1 u2))
	(scale (if (default-object? u2) 1 u1)))
    (assert (and (not (units? value)) (number? scale) (units? units)))
    (if (with-units? value)
	(if (equal? (unit-exponents units)
		    (unit-exponents (u:units value)))
	    value
	    (error "Units do not match: &" value units))
	(with-units (g:* scale (unit-scale units) value)
	  (make-unit (unit-system units)
		     (unit-exponents units)
		     1)))))

So this is actually making a new, wrapped object.

Then, if you look at how it works with the generic system... take addition:

(assign-operation '+   u:+     with-units?             not-differential-or-compound?)
(assign-operation '+   u:+     not-differential-or-compound?  with-units?)

This says that anything that responds true to with-units? can combine with something that responds true to not-differential-or-compound?.

Then we go to definition of u:+:

(define (u:+ x y)
  (cond ((g:zero? x) y)
        ((g:zero? y) x)
        ((units:= x y)
         (with-units (g:+ (u:value x) (u:value y)) (u:units x)))
        ((and *permissive-units*
              (or (without-units? x) (without-units? y)))
         (g:+ (u:value x) (u:value y)))
        (else (error "Units do not match: +" x y))))

And you migt notice that internally, it has a few cases:

• is either side, ie x or y, a zero? if so return the other thing.
• if not, are their units equal? units= if so,
  • unwrap them with (u:value x) etc
  • pass them BACK to generic g:+, where they get re-dispatched and combined
  • re-wrap the result in (with-units result (u:units x))
• the third case is if the variable *permissive-units* is set to true. In that case, then if either side does NOT have units then addition works, but the side with units dumps its units. You might notice that (g:+ (u:value x) (u:value y)) does no rewrapping.
• final case is error, noting that units don't match.

The big pattern with dispatch systems like this is that you can always toss items back out. To simplify a list, for example, (map g/simplify xs) across the list! So good.

How to get dispatch to work in Clojure?

We typically define some type with deftype, then implement the sicmutils.value/Value protocol (see Literal for an example and have it return a
namespaced keyword like ::with-units and ::unit from the kind method.

If you do this, you can write methods like

(defmethod g/add [::v/real ::v/real] [a b]
          (clojure.core/+ a b))

and it will all just work!

[jackrusher]

More references for REPL-optimized symbolic computation with units:

https://reference.wolfram.com/language/tutorial/SymbolicCalculationsWithUnits.html

https://www.gnu.org/software/emacs/manual/html_node/calc/Units.html

[sritchie]

I added a note above about https://painterqubits.github.io/Unitful.jl/stable/

[teodorlu]

Two questions to start off:

1. How do we represent units in data?
2. Are there test cases or examples for the units scheme code?

I'll dig a bit.

[teodorlu]

Scheme examples from units/system.scm:

(with-units->expression SI &foot)
;Value: (& .3048 &meter)

(with-units->expression SI (& 2 &foot))
;Value: (& .6096 &meter)

(with-units->expression SI (/ (* :k (& 300 &kelvin)) :e))
;Value: (& .02585215707677003 &volt)

(with-units->expression SI :c)
;Value: (& 299792458. (* &meter (expt &second -1)))

(with-units->expression SI :h)
;Value: (& 6.6260755e-34 (* (expt &meter 2) &kilogram (expt &second -1)))

[teodorlu]

Should constants constants have a symbolic representation, or should they be normal values? Should we be able to treat c first as a symbol, but be able to substitute the value on demand?

normal:

(require '[sicmutils.constants :refer [c]])

symbolic:

(-> 'c
    (substitute {'c sicmutils.constants/c}))

;; or

(-> 'c
    (substitute (select-keys ['c] sicmutils.constants/defaults)))

[teodorlu]

I'm struggling to find a good starting point. I don't really understand what kind of data model the scheme implementation chose. That makes it hard for me to choose a reasonable first step, because I'm going to be coupled to the data model.

Any suggestions for a starting point for the data model + a reasonable first feature / first unittest?

[sritchie]

@teodorlu , I'll write you tonight, sorry for the delay, lots going on.

I am going to start accumulating some prior art here that we can look at for design, if we want to move beyond / take a different tack than scmutils has taken.

JS

• https://github.com/sharkdp/insect

Clojure

• https://github.com/arrdem/meajure
• https://github.com/mfey/units2
• https://github.com/fogus/minderbinder
• https://github.com/martintrojer/frinj/
• https://github.com/g7s/unit

Julia

This one gets high marks:

• https://painterqubits.github.io/Unitful.jl/stable/

[teodorlu]

Thanks a lot for the writeup! That really helps me with context. I'm looking forward to digging into this. 🤞

[sritchie]

@teodorlu absolutely! Please let me know if you need any help here.

[teodorlu]

@sritchie ,

1. Thanks a bunch for your thorough writeup above. That really helped me get started.
2. I think I'm starting to understand how it all fits together. I'm running your docker image now, evaluating expressions with & and units, which helps unpack how the scheme library works.
3. I'm amazed that (* 'x (& 1 &meter)) is well defined, and that (+ 'x (& 1 &meter)) can error.

I still don't really understand

1. How the dispatch system for the generic operations works; how the "operations with units" can be so seamlessly integrated with the other functions. Perhaps it is as simple as checking for a top-level &, and taking it from there? Not sure.
2. What a Clojure-idiomatic data structure for storing unit information would be. I think a metadata slot might sound reasonable. I'm not sure how that's implemented, whether it's standard Clojure metadata or a sicmutils construct.

Sam's recent introduction of a metadata slot for expressions might be
useful here: could we include a slot in the metadata for units?

We can use the simplifier we have to canonicalize units to their L T M
(etc) form. (in fact our simplifier is overkill because we don't have
addition to worry about.

Not sure what an L T M form is.

[teodorlu]

I'll need to look at sicmutils.generic/def-generic-function dispatching to see how that could work with units.

Ideas:

1. unit info in metadata
2. "high-priority" unit-aware defmethods that handles units, and delegates the math itself to existing generic functions ("dropping" the unit info)

But where does that leave the unit system and quantity (amount + unit) construction? I haven't looked much into the extensibility of scmutils' (Scheme) unit system.

Edit: or do we even want to keep unit info in metadata? Isn't it cleaner to say that a quantity is a value with a unit, and "box in" the value? That avoids dispatch priority conflicts. scmutils (Scheme) seems to take this route, with the & type tag.

[sritchie]

(defn *units [u1 u2]
  (cond (unitless? u1) u2
        (unitless? u2) u1
        :else

        ;; okay, maybe have a nicer error here... NOTE that we can avoid this
        ;; by making our function private: the only way to atually call it
        ;; should be through the dispatch system.
        (do (assert (and (units? u1)
                         (units? u2)))

            ;; make sure they have the same unit system... TODO figure out
            ;; what the unit systems are!
            (assert (and (eq? (unit-system u1)
                              (unit-system u2))))

            (let [v1 (unit-exponents u1)
                  v2 (unit-exponents u2)]
              (cond (empty? v1)
                    (make-unit (unit-system u1)
                               v2
                               (* (unit-scale u1)
                                  (unit-scale u2)))

                    (empty? v2)
                    (make-unit (unit-system u1)
                               v1
                               (* (unit-scale u1)
                                  (unit-scale u2)))

                    :else (make-unit
                           (unit-system u1)
                           (mapv + v1 v2)
                           (* (unit-scale u1)
                              (unit-scale u2))))))))

[sritchie]

(deftype Unit [exponents v scale] ,,)

[teodorlu]

working title: should we support late-bound units?

Scmutils treats units for numbers and symbols the same: either all units are set in advance, or there's no units. It cannot await the definition of the unit for a symbol.

This surprised me, because I was expecting to substitute a number and a unit for a symbol.

A number can be added to a symbol, provided they are of the same unit (in advance):

1 ]=> (+ (& 1 &meter) (& 'x &meter))
#|
(& (+ 1 x) &meter)
|#

If we avoid setting a unit on x, and add to 1 meter, we crash:

1 ]=> (+ (& 1 &meter) 'x)

;Units do not match: +
;To continue, call RESTART with an option number:
; (RESTART 1) => Return to read-eval-print level 1.

Should we allow to substitute the expression x with x=1 m, or x m with x=1?

---

1. This is not a blocker for sicmutils/sicmutils#213, and is not urgent. Hammock time candidate.
2. For making a choice about better support for "more flexible SI units", I think it would be good to have a motivating example from dimensional analysis. What is a specific example where the scheme implementation is not sufficient?

---

ref

https://github.com/sicmutils/sicmutils/blob/cb527804c274448d75a2d467078c89246c6aa0aa/src/sicmutils/units/constants.cljc#L5-L31

[teodorlu]

Motivation for late bound units: support Dimensional analysis

[teodorlu]

How should units be printed?

In Scheme, there's global namespacing. scmutils simply prints meters as &meter, where &meter is available globally.

That doesn't work with Clojure namespaces.

EDN reader tags / data readers could be a solution. The time-literals library does this for time types: https://github.com/henryw374/time-literals/blob/master/src/data_readers.cljc

data_readers.cljc, when found on the classpath, seems to make reader tags globally available. But data_readers.cljc is static. So it could provide sicmutils/si-unit.

But if we want scheme style custom unit systems, we also need to provide tags like this dynamically. If you create your-ns/your-custom-si-system, then we'd need to generate a reader literal for that. Which probably means def-unit-system needs to be a macro.

[teodorlu]

Custom unit systems does sound a little scary. I'm not that comfortable writing macros.

Starting out with providing working SI units sounds like a good idea to me. Though let's not assume there's ever only going to be one unit system.

[teodorlu]

Can we multiply numbers from different unit system?

scmutils says no:

1 ]=> (define-unit-system 'teod (list '&dollar "$" "dollar"))
#| teod |#
1 ]=> (* 1 &dollar)
#| &dollar |#
1 ]=> (* 1 &dollar &meter)
;Assertion failed: (and (eq? (unit-system u1) (unit-system u2)))
;To continue, call RESTART with an option number:
; (RESTART 1) => Return to read-eval-print level 1.

via Slack: https://clojurians.slack.com/archives/C01ECA9AA74/p1621075975042800

[jackrusher]

How's this work going, lads?

[teodorlu]

Heh.

I figured that I had gotten in over my head. This is harder than other problems I've solved with Clojure. I haven't put any serious effort for SI Units in Sicmutils this year.

If I were to try again, I'd isolate the unit problem. Try to figure out how to provide "just units in a nice way" with plain Clojure data. Find a nice API for units. Then tackle integration with sicmutils.

Two things I couldn't get my head around the first time:

1. Should we hard-code SI Units as the only unit system? If not, how do we handle different unit systems?
2. I'd like this to be nice to use from a REPL. I think a reader literal could be a nice solution. Something like #with-unit [30 kN m].

[teodorlu]

@2food released a Clojure library for SI units today: https://github.com/anteoas/broch