- home
- ||
- publications
- ||
- blog
- ||
- contact

Dimensional analysis is quite the tool to track down stupid
mistakes in physics much like a type-checker will detect non-sensical
expressions in your favourite statically-typed programming language.
I have been meaning to toy around with this notion in a dependently-typed
setting for quite a while. Here are a few definitions. As well as an
Agda program that needs to be *compiled*.

On Wednesday, Bob presented at SPLS a really cool project he has been working on: a type system which, by parametricity, derives automagically [1] the premises of Noether's theorem. All this talking about physics got me motivated enough to spend Saturday's afternoon hacking up a toy example. This is probably not revolutionary [2] but it was quite amusing to develop anyway so here it is.

This whole project is structured as a succession of modules in order to to have separate name spaces. This way, we can make sure that the different operations morally implementing the same concepts can have the same name.

We will limit ourselves to 3 types of units of measure here: kilograms, meters and seconds. The dimension of an object is modelled by a record storing the exponents assigned to each one of these components.

record dimension : Set where field kilogram : ℤ meter : ℤ second : ℤ

Taking the product of two dimensions amounts to summing the degrees for each one of the units of measure. Quotienting is, likewise, a pointwise operation on the two vectors: this time we compute the difference rather than the sum.

_*_ : (d e : dimension) → dimension d * e = record { kilogram = kilogram d + kilogram e ; meter = meter d + meter e ; second = second d + second e }

We can now define the basic dimensions: kilograms, meters and seconds by assigning + 0 to all the fields except for the one of interest which is instantiated with + 1. As an example, here is the dimension sec corresponding to seconds.

sec = record { kilogram = + 0 ; meter = + 0 ; second = + 1 }

So far, so good. But in real life we quite like to consider other units than just the basic ones. A unit of measure is therefore a dimension together with a non-zero coefficient. The type of hn has been chosen so hn can be inferred by the system in concrete cases (if n is non-zero, then hn is of type ⊤ hence equal to tt) and to match the requirements imposed by _div_ in the standard library.

data unit : Set where _,_#_ : (n : ℕ) (hn : False (n ℕ.≟ 0)) (d : dimension) → unit

The product of two units is once more a pointwise operations: it is the combination of the product of their coefficients (which is guaranteed to be non-zero) and the product of their respective dimensions.

_*_ : (u v : unit) → unit (k , hk # d) * (l , hl # e) = k * l , _ # d * e

It is quite natural to introduce the SI's prefixes as functions altering a previously defined unit. In our system, the types of milli, centi and deci are nastier than the ones presented below because we are using natural numbers rather than rational ones and division of a non-zero number by a non-zero number can return zero.

deca hecto kilo : unit → unit

A minute is nothing more than 60 seconds. 60 being different from 0 (!), the proof obligation reduces to the only inhabitant of unit which is inferred by the system:

min : unit min = 60 , _ # sec

Finally, values are just numerical values repackaged together with a unit. We choose here to have a data constructor forcing the user to explicitly mention the unit of measure to be attached to the value. It could be avoided given that d is a parameter of the data-type but it is quite handy for documentation purposes and on the fly definition of values.

data \[_] : (d : unit) → Set where ⟨_∶_⟩ : (value : ℕ) (d : unit) → \\[ d ]

Multiplying to value has to return something of the right dimension but the scaling factor can be arbitrary (we may multiply kilo meters by hours and expect a result in m / s). Hence the following implementation. Addition and division behave similarly.

_*_ : \\[ k , hk # d ] → \\[ l , hl # e ] → \\[ m , hm # d * e ] _*_ ⟨ vd ∶ ._ ⟩ ⟨ ve ∶ ._ ⟩ = ⟨ (k * vd * l * ve) div m ∶ _ ⟩

Our implementation of value multiplication is so generic in its type that it does not have the subformula property which will probably be problematic in large expression where the type of all the subexpressions is not inferrable. A simple solution is to specialize the return type down to \\[ 1 , _ # d * e ] whilst introducing a lifting operation converting from one type to an other.

↑ : ∀ {k hk l hl d} → \[ k , hk # d ] → \[ l , hl # d ] ↑ ⟨ v ∶ ._ ⟩ = ⟨ (k * v) div l ∶ _ ⟩

This system tracking dimensions makes sure that we do not make mistakes when combining different values but it also converts between various units. For instance, the following expression normalizes to 10.

60hm/min : \[ deca m / s ] 60hm/min = ⟨ 60 ∶ hecto m ⟩ / ⟨ 1 ∶ min ⟩

The setting of this simulation is pretty simple: we have a ball
characterized by its position, speed and acceleration in a vertical plan.
The gravitational constant g corresponding to earth's attraction
is expressed in m / (s * s) and equal to 10.

In order to simulate the free fall of this ball, we are going to apply
Newton's approximation method: time is discretized and the characterization
of the ball is updated at every single time-step. Here the type system
guarantees that we do not make mistakes.

newton : ∀ (dt : \[ s ]) (p : point) → point newton dt p = record { accx = accx p ; accy = accy p + g ; vx = vx p + (accx p :* dt) ; vy = vy p + (accy p :* dt) ; x = x p + (vx p :* dt) ; y = y p + (vy p :* dt) }

A simulation consists in applying repeatedly this newton transformation to an initial point. We define the throw function doing exactly this and generating a trace of all the successive positions.

throw : (n : ℕ) (dt : \[ s ]) → point → Vec point (ℕ.suc n) throw zero dt p = p ∷ \[] throw (suc n) dt p = p ∷ throw n dt (newton dt p)

After a little bit of magic, we are able to call a tracing function from Haskell's Gnuplot wrapper and generate the following graph showing the (upside-down) ball's trajectory in the vertical plan.

First of all, it works quite well for a first experiment! In a serious
implementation though, one should really use, at the very least, rational
numbers to represent the coefficients. This would ensure that conversion
ratios can be properly computed and make the work a lot more useful.
It is no coincidence that my example involving a speed expressed in
*distance per minutes* precisely picks a distance which is a multiple of
60... I haven't had the internet connection, the time nor the
desire to see whether appropriate libraries already exist for Agda or to
write my own ones.

Similarly, the fact that the free fall is represented upside-down is a
direct consequence of the values not being integers but rather natural
numbers (once more due to a lack of libraries in Agda: div and
mod are only available for ℕ in the standard library).

And Agda's automatic case-splitting, support for Unicode, definition of
functions *via* equations, my (ridiculously small but good enough)
understanding of how to import Haskell function as postulates, etc. kept me
away from Coq (at least for the experimentation phase).

A good start to learn about the free lunch
provided by parametricity is obviously Phil Wadler's *Theorems for Free!*
(citeulike).

Andrew Kennedy's paper *Types for
Units-of-Measure: Theory and Practice* mentions that abusing C++ or
Haskell's type systems makes it possible to « achieve it, but at some cost
in usability ». Arguably, a real dependently-typed language is more suited
for computing at the type level in this fashion.

Last update: 2016 11 24