Scientific.Dimension: Type Arithmetic and Physical Units in Haskell

June 3, 2007

One of Haskell’s most outstanding features is it’s type-checker. Mis-matching types are a quick indication that something is wrong with a particular piece of code. And type annotations are useful documentation tool that provides hints about the uses of a function. Type checking is however not unique to programming languages. There is a similar concept in Physics: Units. Some programming languages, such as Mathematica provide built-in unit-checking. The Haskell language lends itself to this purpose very nicely, due to it’s type-checking system. As opposed to Mathematica, this can be done without changes to the core language itself.

This post will start out by introducing simple type arithmetic using Peano numerals and progresses to implementing a full unit checker on top of the Haskell language. The post is also available as a literate Haskell module.

The Glasgow Haskell Compiler (GHC) is used, as we will require the use of the -fglasgow-exts and -fallow-undecidable-instances extensions.

> {-# OPTIONS -fglasgow-exts #-}
> {-# OPTIONS -fallow-undecidable-instances #-}
> module Scientific.Dimensions where

Arithmetic in the type-checker

In our unit-checking code, we will need the ability to add and subtract numbers in the type system. The numbers themselves are represented using Peano numerals. Positive Peano numerals are composed of two parts: Succ and Zero. In Haskell, their definition would look like this:

> -- Peano numerals
> data Zero = Zero deriving (Eq, Show)
> data Succ n = Succ n deriving (Eq, Show)

Where Zero represents the number “0”, and Succ is the “successor” of a number. For example, the type ‘Succ Zero’ is “1”, and ‘Succ (Succ (Succ Zero))’ is the number “3”. To add two numbers, we create a Haskell type class.

> class Add a b c | a b -> c where
>   add :: a -> b -> c

In the above class definition, the type variables ‘a’ and ‘b’ are the inputs to the calculation, and ‘c’ is the resulting output type. The “a b -> c” part tells the type-checker that the type-variable ‘c’ can be determined by knowing ‘a’ and ‘b’. In other words, a and b “determine” c. This means that for every possible case of ‘a’ and ‘b’, we will have to specify the resulting ‘c’. Here’s our first case:

> -- Adding any number to Zero returns that number
> instance Add num Zero num where
>   add num Zero = num

Then, we need to cover the reverse condition: adding Zero to a positive number. We only accept positive numbers here, since the case of adding Zero to Zero is already covered by the previous case.

> -- Adding Zero to a positive number returns that number
> instance Add Zero (Succ a) (Succ a) where
>   add Zero (Succ a) = Succ a

Now, we need to cover the slightly more complicated case of adding any two positive numbers. This is done recursively: Given the positive numbers ‘Succ a’ and ‘Succ b’, we first calculate ‘c = a + b’, and then return ‘c + 2′. Note that ‘a’ is equal to ‘Succ a’ minus one, and that ‘c + 2′ is the same as ‘Succ (Succ c)’.

> -- Adding two positive numbers 'Succ a' and 'Succ b' returns 'Succ (Succ (a + b))'
> instance Add a b c => -- This means: c = a + b
>          Add (Succ a) (Succ b) (Succ (Succ c)) where
>   add (Succ a) (Succ b) = Succ (Succ c)
>                 where c = add a b

In the above instance declaration, everything before the ‘=>’ symbol is a precondition. Since we calculate ‘c = add a b’, we need to make sure that addition is defined on ‘a’ and ‘b’ with result-type ‘c’. This is the reason for the pre-condition.

In our unit-checking code, we will also need negative numbers. They are defined using ‘Prec’, where ‘Prec’ stands for “preceding”. For example, ‘Prec (Prec Zero)’ represents the number “-2″.

> -- Negative Peano numerals
> data Prec n = Prec n

Of course, we will have to extend our definition of ‘Add’ to include negative numbers. The new instance declarations are very similar to the ones for positive numbers.

> -- adding Zero to a negative number returns that number
> instance Add Zero (Prec a) (Prec a) where
>   add Zero (Prec a) = Prec a
> -- Adding negative numbers 'Prec a' and 'Prec b' returns 'Prec (Prec (a + b))' 
> instance Add a b c =>
>          Add (Prec a) (Prec b) (Prec (Prec c)) where
>   add (Prec a) (Prec b) = Prec (Prec c)
>                 where c = add a b
> -- Adding negative number 'Prec a' to positive number 'Succ b' returns 'a + b'
> instance Add a b c =>
>          Add (Prec a) (Succ b) c where
>   add (Prec a) (Succ b) = add a b
> -- Reverse of above
> instance Add a b c =>
>          Add (Succ a) (Prec b) c where
>   add (Succ a) (Prec b) = add a b

That’s it! We are now able to add any two numbers. Let’s try it out:

> -- 1 = 3 - 2
> one = add (Succ (Succ (Succ Zero))) (Prec (Prec Zero))

The type of ‘one’ is ‘Succ Zero’. This is determined at compile-time using our recursive class-definition for ‘Add’. We will also want to be able to subtract two numbers. The definition of Sub is almost identical to Add:

> class Sub a b c | a b -> c where
>   sub :: a -> b -> c
> instance Sub a Zero a where
>   sub a Zero = a
> instance Sub Zero (Succ a) (Prec a) where
>   sub Zero (Succ a) = Prec a
> instance Sub Zero (Prec a) (Succ a) where
>   sub Zero (Prec a) = Succ a
> instance Sub a b c =>
>          Sub (Succ a) (Prec b) (Succ (Succ c)) where
>   sub (Succ a) (Prec b) = Succ (Succ c)
>                 where c = sub a b
> instance Sub a b c =>
>          Sub (Prec a) (Succ b) (Prec (Prec c)) where
>   sub (Prec a) (Succ b) = Prec (Prec c)
>                 where c = sub a b
> instance Sub a b c =>
>          Sub (Prec a) (Prec b) c where
>   sub (Prec a) (Prec b) = sub a b
> instance Sub a b c =>
>          Sub (Succ a) (Succ b) c where
>   sub (Succ a) (Succ b) = sub a b

Lists in the Type System

We will also need to use lists in our type-checker. Lists have a very simple definition:

> -- Linked list
> data Nil = Nil deriving (Eq, Show)
> data x `Cons` xs = x `Cons` xs deriving (Eq, Show)

Nil is an empty list, and Cons is a pair of types. A list of an Integer, a String, and a Double is defined like this:

> funkyList :: Int `Cons` (String `Cons` (Double `Cons` Nil))
> funkyList = 256 `Cons` ("hello world" `Cons` (3.14159 `Cons` Nil))

Changing the fixity of the Cons type-constructor will allow us to drop the parentheses.

> infixr 6 `Cons`

Using this definition of lists allows many functions to be implemented in a more type-safe manner than using Haskell lists. For example, the Prelude functions ‘head’ and ‘tail’ will fail when passed an empty list. If head was instead defined using our new definition for lists, the type-checker would ensure at compile-time that no empty lists are passed as arguments.

> safeHead :: x `Cons` xs -> x
> safeHead (x `Cons` xs) = x

Running something like ‘safeHead Nil’ will result in a type-error. Implementing functions such as ‘foldl’, ‘map’ and ‘concat’ is left as an exercise to the reader. (If you can’t figure it out, make sure to ask for help at #haskell on irc.freenode.org)

Unit Checking

We are now ready to start start implementing our unit-checker. Units are represented as a list of numbers that count how much of a certain “primitive” unit is contained in our unit. For example, suppose we use use the SI units Meter, Kilogram and Second as our primitive units:

> -- primitive unit count
> data MeterCount = MeterCount deriving (Eq, Show)
> data KilogramCount = KilogramCount deriving (Eq, Show)
> data SecondCount = SecondCount deriving (Eq, Show)

The unit Meter would be defined as having a MeterCount of “1”, a KilogramCount of “0” and a SecondCount of “0”:

> type Meter = (MeterCount, Succ Zero)
>       `Cons` (KilogramCount, Zero)
>       `Cons` (SecondCount, Zero)
>       `Cons` Nil

Recall that ‘Succ Zero’ is the Peano numeral representing “1”. Similarly, the unit of ‘meter^2 * kilogram / second^3′ would be defined as:

> type FunkyUnit = (MeterCount, Succ (Succ Zero))
>           `Cons` (KilogramCount, Succ Zero)
>           `Cons` (SecondCount, Prec (Prec Zero))
>           `Cons` Nil

Kilograms and Seconds are defined similar to Meters:

> type Kilogram = (MeterCount, Zero)
>          `Cons` (KilogramCount, Succ Zero)
>          `Cons` (SecondCount, Zero)
>          `Cons` Nil
> type Second = (MeterCount, Zero)
>        `Cons` (KilogramCount, Zero)
>        `Cons` (SecondCount, Succ Zero)
>        `Cons` Nil

Unitless values occur quite frequently in physics, so we create a type for them as well:

> type Unitless = (MeterCount, Zero)
>          `Cons` (KilogramCount, Zero)
>          `Cons` (SecondCount, Zero)
>          `Cons` Nil

To multiply two units, we define the Mult class:

> -- multiplies to units
> class Mult a b c | a b -> c, a c -> b, b c -> a where
>   mult :: a -> b -> c

As with the Add class, ‘a’ and ‘b’ are the two units multiplied together, and ‘c’ is the resulting unit. This time, any two units determine the third. As with the ‘Add’ class, we need to instantiate the ‘Mult’ class with every possible input.

> -- multiply empty units
> instance Mult Nil Nil Nil where
>   mult Nil Nil = Nil
> -- multiply units composed of non-empty lists recursively
> instance (Add i j k, Mult r s t) =>
>          Mult ((a, i) `Cons` r) ((a, j) `Cons` s) ((a, k) `Cons` t) where
>   mult ((a, i) `Cons` r) ((_, j) `Cons` s) h1>
>     (a, add i j) `Cons` mult r s

Division of units is very similar to multiplication. In stead of adding unit-counts, they are subtracted.

> class Div a b c | a b -> c, a c -> b, b c -> a where
>   divide :: a -> b -> c
> -- divide empty units
> instance Div Nil Nil Nil where
>   divide Nil Nil = Nil
> -- divide units composed of non-empty lists recursively
> instance (Sub i j k, Div r s t) =>
>          Div ((a, i) `Cons` r) ((a, j) `Cons` s) ((a, k) `Cons` t) where
>   divide ((a, i) `Cons` r) ((_, j) `Cons` s) h1>
>     (a, sub i j) `Cons` divide r s

Now, we create a Dimension type: a value-unit pair.

> -- Physical dimension with unit
> data value `Dimension` unit = value `Dimension` unit deriving (Eq, Show)

We also want `Dimension` to have lower precedence than `Cons` to eliminate unnecessary parentheses:

> infix 4 `Dimension`

To make it easy to create physical dimensions, we create constructors for each of our primitive unit types.

> meter :: Fractional x => x -> x `Dimension` Meter
> meter x = x `Dimension` (MeterCount, Succ Zero)
>	           `Cons` (KilogramCount, Zero)
>		   `Cons` (SecondCount, Zero)
>      		   `Cons` Nil
> gram :: Fractional x => x -> x `Dimension` Kilogram
> gram x = (0.001 * x) `Dimension` (MeterCount, Zero)
>      	               `Cons` (KilogramCount, Succ Zero)
>      	 	       `Cons` (SecondCount, Zero)
>      		       `Cons` Nil
> second :: Fractional x => x -> x `Dimension` Second
> second x = x `Dimension` (MeterCount, Zero)
>      	            `Cons` (KilogramCount, Zero)
>      	  	    `Cons` (SecondCount, Succ Zero)
>      		    `Cons` Nil
> unitless :: Fractional x => x -> x `Dimension` Unitless
> unitless x = x `Dimension` (MeterCount, Zero)
>      	              `Cons` (KilogramCount, Zero)
>      		      `Cons` (SecondCount, Zero)
>      		      `Cons` Nil

For SI units it makes sense to define unit-prefixes (kilo, cent, milli, etc.):

> giga f x  = f (x * 1000000000)
> mega f x  = f (x * 1000000)
> kilo f x  = f (x * 1000)
> deca f x  = f (x * 10)
> deci f x  = f (x * 0.1)
> centi f x = f (x * 0.01)
> milli f x = f (x * 0.001)
> micro f x = f (x * 0.000001)
> nano f x  = f (x * 0.000000001)

Using the primitive dimension constructors and prefixes, dimensions are defined like thus:

> mile = meter 1609.3
> gravity = kilo gram 9.8 .* meter 1 ./ (second 1 .* second 1)
> speedOfLight = unitless (10^8) .* meter 3 ./ second 1

It would be nice, if we could instantiate the Prelude Num class to define functions such as (+), (-) and (*) for our dimensions. However, multiplication in the Num class is only allowed on values of equal type. Thus, we need to define our own addition, subtraction, multiplication and division functions:

> (a `Dimension` b) .* (c `Dimension` d) = (a * c) `Dimension` (mult b d)
> (a `Dimension` b) ./ (c `Dimension` d) = (a / c) `Dimension` (divide b d)

For addition and subtraction, we need to additionally make sure that the units of the two arguments are the same:

> (.+) :: Num a => (a `Dimension` b) -> (a `Dimension` b) -> (a `Dimension` b)
> (a `Dimension` b) .+ (c `Dimension` _) = (a + c) `Dimension` b
> (.-) :: Num a => (a `Dimension` b) -> (a `Dimension` b) -> (a `Dimension` b)
> (a `Dimension` b) .- (c `Dimension` _) = (a - c) `Dimension` b

However, it does make sense for Unitless dimensions to instantiate Num. That’s exactly what unitless dimensions are: numbers without unit.

> instance Num a => Num (a `Dimension` Unitless) where
>   (a `Dimension` unit) + (b `Dimension` _) = (a + b) `Dimension` unit
>   (a `Dimension` unit) - (b `Dimension` _) = (a - b) `Dimension` unit
>   (a `Dimension` unit) * (b `Dimension` _) = (a * b) `Dimension` unit
>   abs (a `Dimension` unit) = abs a `Dimension` unit
>   signum (a `Dimension` unit) = signum a `Dimension` unit
>   fromInteger i = fromInteger i `Dimension` unit
>     where _ `Dimension` unit = unitless 1

To use the (/) function on Unitless values, we need to instantiate the ‘Fractional’ class:

> instance Fractional a => Fractional (a `Dimension` Unitless) where
>   (a `Dimension` unit) / (b `Dimension` _) = (a / b) `Dimension` unit
>   fromRational r = fromRational r `Dimension` unit
>     where _ `Dimension` unit = unitless 1

Conclusion

Unit checking is very useful for Scientific applications, and is now possible in the Haskell language. Hopefully, this will encourage Haskell’s use in the Scientific community.

About these ads

6 Responses to “Scientific.Dimension: Type Arithmetic and Physical Units in Haskell”


  1. […] Scientific.Dimension: Type Arithmetic and Physical Units in Haskell One of Haskell’s most outstanding features is it’s type-checker. Mis-matching types are a quick indication […] […]

  2. Ossi Herrala Says:

    I think you are missing couple of zeroes here:

    > giga f x = f (x * 1000) — x * 1000000000
    > mega f x = f (x * 1000) — x * 1000000
    > kilo f x = f (x * 1000)

  3. liftm Says:

    Ossi: Thanks for the correction.


  4. I’m aiming to provide a reasonably complete library with statically unit checking. The project site is http://code.google.com/p/dimensional/. I’m currently trying to wrap the library up for a 1.0 release.

    It’s good to see that someone else is interested in this “problem”. While there are some significant differences (e.g. I’ve differentiated between units and “quantities” at the type level) I think we’ve tackled it similarly.

    Feel free to take a look at my library. The main documentation is the extensively commented literate haskell code and given that you have solved the problem once it should be easy for you to follow. I imagine e.g. the support for powers and roots might interest you.

    If your interest in this goes beyond a one-off academic exercise maybe you would be interested in contributing/collaborating? Thanks.

  5. Mietek Bąk Says:

    There’s a problem with your Sub instances, as evidenced by the SecondCount of your gravity example being Prec (Succ Zero).

    This should fix the problem:

    > instance (Sub Zero a b) => Sub Zero (Succ a) (Prec b) where
    > sub Zero (Succ a) = Prec (sub Zero a)
    >
    > instance (Sub Zero a b) => Sub Zero (Prec a) (Succ b) where
    > sub Zero (Prec a) = Succ (sub Zero a)

  6. sandrar Says:

    Hi! I was surfing and found your blog post… nice! I love your blog. :) Cheers! Sandra. R.


Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s

Follow

Get every new post delivered to your Inbox.

%d bloggers like this: