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.

# The data structures involved

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.

# Dimensions

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 }

# Units of measure

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

# Values

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 ∶ _ ⟩

# Examples

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 ⟩

# Application: simulating the free fall of a ball

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.

# Conclusions

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

# Footnotes