Compile time dimensional analysis with Libra

Dimensional analysis

When we code, we code in numbers - doubles, floats and ints. Those numbers always represent real world quantities.

For example, the number of people in a room can be represented as an integer, as can the number of chairs. Adding people and chairs together gives a nonsensical result, but dividing the number of people by the number of chairs gives a useful indicator of how full up the room is.

val numberOfPeople = 9
val numberOfChairs = 10
numberOfPeople + numberOfChairs // this is a bug
numberOfPeople.toDouble / numberOfChairs.toDouble // this is useful

This is actually a form of dimensional analysis. We're mentally assigning the dimension Person to the quantity of people, and Chair to the quantity of chairs. Dimensional analysis can be summarized in two laws.

1. Quantities can only be added or subtracted to quantities of the same dimension
2. Quantities of different dimensions can be multiplied or divided

Why is it important?

Ignoring the laws can result in serious problems. Take the Mars Climate Orbiter, a $200 million space probe which successfully reached Mars after a year long voyage, but suddenly crashed into the Martian atmosphere on arrival. Most components on the orbiter were using metric units, however a single component was sending instructions in Imperial units. The other components did not detect this, and instead began a sudden descent causing the orbiter to burn up. This was a simple unit conversion error! It was a basic mistake that could have been easily avoided. It should have been picked up during testing, or in the runtime validation layer.

In fact, it could even have been caught at compile time.

Compile time dimensional analysis

We're going to use a similar problem to demonstrate compile time dimensional analysis. To fit with the theme of rocket physics, we will tackle a rocket launch towards the distant constellation of Libra. We'll begin by working through our calculation in doubles before adding compile time safety with dependent types and finally supporting compile time dimensional analysis with typeclass induction.

Destination: Alpha Librae

The star that we're aiming for is Alpha Librae. This is pretty far, so we can only send one very small person. We have been given the following quantities to work with:

• rocket mass of a small person - $40 \text{kg}$
• fuel mass of a lot of fuel - $10^4 \text{kg}$
• exhaust speed of a decent fuel - $10^6 \text{ms}^{-1}$
• distance to Alpha Librae - $77 \text{ly}$

We want to calculate when the rocket will arrive.

To do so, we're going to make use of a formula known as the Ideal Rocket Equation. This calculates the speed of a rocket in ideal conditions.

val rocketSpeed = exhaustSpeed * log((rocketMass + fuelMass) / rocketMass)

Once we have the speed, we can work out the travel time.

val time = distance / rocketSpeed

Plugging the numbers in

Let's do what we're used to doing and use doubles:

val rocketSpeed = 1000000.0 * log((40.0 + 10000.0) / 40.0)
// rocketSpeed: Double = 5525452.939131783

val time = 77.0 / rocketSpeed
// time: Double = 1.39355091515989E-5

Fantastic! We can get to Libra in less than a day!

Unfortunately, this time estimate is too far off to be valid. We can't get to Libra that quickly at light speed, let alone rocket speed. We've clearly made a mistake somewhere. Instead of pouring over our code to find out where that is, let's try and use the compiler.

Using types

We can add some type safety to this problem by using a case class to represent each quantity.

case class Quantity[A](value: Double)

A represents the quantity dimension. So given the following dimensions:

type Kilogram
type Metre
type Second
type MetresPerSecond
type C
type LightYear
type Year

We can create quantities:

val rocketMass = Quantity[Kilogram](40.0)
val fuelMass = Quantity[Kilogram](10000.0)
val exhaust = Quantity[MetresPerSecond](1000000.0)
val distance = Quantity[LightYear](77.0)

It's important to note that these are types, not classes. We never instantiate a MetresPerSecond - we're just using it to differentiate between Quantity[MetresPerSecond] and Quantity[Year] at the type level.

So how does this change the code?

val rocketSpeed = Quantity[MetresPerSecond](exhaust.value * log((rocketMass.value + fuelMass.value) / rocketMass.value))
// rocketSpeed: Quantity[Types.MetresPerSecond] = Quantity(5525452.939131783)

val time = Quantity[Year](distance.value / rocketSpeed.value)
// time: Quantity[Types.Year] = Quantity(1.39355091515989E-5)

In short, it doesn't. The code might be clearer, but we don't know what the bug is. This is because the compiler isn't doing anything with the types we've added.

Operating on quantities

We can encode our first law of addition at compile time by creating a function to add quantities:

def add[A](q0: Quantity[A], q1: Quantity[A]): Quantity[A] = Quantity(q0.value + q1.value)

This ensures that quantities can only be added to other quantities of the same type. Trying to add quantities of different types will result in a compilation error.

A quantity can also be multiplied by a dimensionless scalar value to give a quantity of the same dimension.

def times[A](q: Quantity[A], v: Double): Quantity[A] = Quantity(q.value * v)

It would be great if we could divide quantities too. Writing a divide function is more difficult:

def divide[A, B, Magic](q0: Quantity[A], q1: Quantity[B]): Quantity[Magic] =
  Quantity[Magic](q0.value / q1.value)

There's a clear problem with trying to do this. When we divide a quantity by another, we don't know what the Magic output type should be. The output type is dependent on what the input types are (for example, dividing Metre by Second should give MetresPerSecond). The compiler needs a way of working out what the output is, provided that it knows the input types.

Dependent types

What we actually want is a dependent type. A division operation should occur at the type level, taking two input types and supplying a dependent output type. We can create the trait Divide with a dependent output type:

trait Divide[A, B] {
  type Out
}

We also need to define an Aux type alias. This is known as the Aux pattern and makes it easier to refer to all three types at once.

object Divide {
  type Aux[A, B, Out0] = Divide[A, B] { type Out = Out0 }
}

We can create instances of this divide typeclass with different output types, so the output type is dependent on the value of the divide typeclass instance.

When dividing, the compiler looks for this implicit typeclass instance and returns a quantity corresponding to the output type.

def divide[A, B](q0: Quantity[A], q1: Quantity[B])(implicit d: Divide[A, B]): Quantity[d.Out] =
  Quantity[d.Out](q0.value / q1.value)

So given that we want to divide A by B, the compiler will look for a value of Divide[A, B] and find the Out type of it. If no instance exists, the code doesn't compile.

We'll need some more types to represent the result of a division:

type LightYearSecondsPerMetre
type MetresPerSecondPerC
type Dimensionless

And we'll need to write instances for all combinations of dimensions.

implicit val kgDivideKg: Divide.Aux[Kilogram, Kilogram, Dimensionless] =
  new Divide[Kilogram, Kilogram] { type Out = Dimensionless }

implicit val lyDivideC: Divide.Aux[LightYear, C, Year] =
  new Divide[LightYear, C] { type Out = Year }

implicit val lyDivideMps: Divide.Aux[LightYear, MetresPerSecond, LightYearSecondsPerMetre] =
  new Divide[LightYear, MetresPerSecond] { type Out = LightYearSecondsPerMetre }

implicit val mpsDivideC: Divide.Aux[MetresPerSecond, C, MetresPerSecondPerC] =
  new Divide[MetresPerSecond, C] { type Out = MetresPerSecondPerC }

implicit val mpsDivideMpsPerC: Divide.Aux[MetresPerSecond, MetresPerSecondPerC, C] =
  new Divide[MetresPerSecond, MetresPerSecondPerC] { type Out = C }

And so on.

Unfortunately, there are an infinite number of combinations, so there are an infinite number of instances. Nevertheless, let's plough on with the ones we've written. We can modify our rocket equation to use add, times and divide:

val rocketSpeed = times(exhaust, log(divide(add(rocketMass, fuelMass), rocketMass).value))
// rocketSpeed: Quantity[Types.MetresPerSecond] = Quantity(5525452.939131783)

val time: Quantity[Year] = divide(distance, rocketSpeed)
// <console>:31: error: type mismatch;
//  found   : Quantity[lyDivideMps.Out]
//     (which expands to)  Quantity[MoreTypes.LightYearSecondsPerMetre]
//  required: Quantity[Types.Year]
//        val time: Quantity[Year] = divide(distance, rocketSpeed)
//                                         ^

Great! We've caught our bug! The result was in LightYearSecondsPerMetre, not Year. We made a unit conversion error, just like the Mars orbiter.

We can now fix this by adding a conversion:

val metresPerSecondPerC: Quantity[MetresPerSecondPerC] = divide(Quantity[MetresPerSecond](300000000.0), Quantity[C](1.0))
// metresPerSecondPerC: Quantity[MoreTypes.MetresPerSecondPerC] = Quantity(3.0E8)

val speedInC = divide(rocketSpeed, metresPerSecondPerC)
// speedInC: Quantity[mpsDivideMpsPerC.Out] = Quantity(0.018418176463772612)

val time: Quantity[Year] = divide(distance, speedInC)
// time: Quantity[Types.Year] = Quantity(4180.65274547967)

It seems like it's going to take a lot longer than we hoped to get to Libra. Perhaps it's unwise to send a person.

Automatic derivation

We found the bug, but we needed to explicitly write out typeclass instances for every combination of dimensions. This might have worked for our small problem, but it just doesn't scale in the long run. We need to figure out a way of deriving the typeclass instances automatically. To attempt this, we first need to generalize what a combination of dimensions actually is.

Representing dimensions

We can represent a combination of dimensions as a heterogeneous list (HList) of base dimensions. HLists are defined in shapeless, a cornerstone of most functional libraries, and can be thought of as a type level list.

type LightYearSeconds = LightYear :: Second :: HNil

This is good for multiples of dimensions, such as LightYearSeconds, but doesn't represent combinations created from division, such as MetresPerSecond. To do this, we need some way of representing integer exponents as types. We can represent integers as types using Singleton types. We actually need these singleton types in type position. This is supported by a new feature present in Typelevel Scala, called literal types:

scalaOrganization := "org.typelevel"
scalacOptions += "-Yliteral-types"

We need to represent a key value pair of dimension and integer exponent. We could use a Tuple for this, but will use a shapeless FieldType instead. This is similar to a Tuple, but is more compatible with some of shapeless's typeclasses.

type MetresPerSecond = FieldType[Metre, 1] :: FieldType[Second, -1] :: HNil

It's important to note that the number 1 above is a type, not a value. Because it's a type, the compiler can work with it.

Operations on Singleton types

When we multiply and divide dimensions, we want to add or subtract from these exponents.

We can use a library called singleton ops to do this. This provides us with type level integer operations using the OpInt typeclass:

OpInt.Aux[1 + 2, 3]
OpInt.Aux[3 * 2, 6]

It also provides a convenient alias for integer singleton types

type XInt = Singleton with Int

The type 1, for example, is a subtype of XInt.

Deriving typeclass instances

We now need to automatically derive typeclass instances of Divide. To do this, we're going to derive instances for Invert and Multiply operations first. Deriving Divide then becomes much simpler.

The technique we're going to use to automatically derive instances is known as typeclass induction.

Typeclass Induction

Aaron Levin gave a great introduction to induction in his talk earlier at the Typelevel Summit. In summary, you can derive an implicit typeclass instance for all cases by:

1. Providing it for the base case
2. Providing it for the n + 1 case, given that the n case is provided

This is similar to the mathematical method of proof by induction.

Invert

We're first going to derive inductive typeclass instances for the Invert operation. Inverting a quantity raises it to the exponent of -1. This means that the exponents of all dimensions must be negated.

For example, the inverse of FieldType[Metre, 1] :: HNil is FieldType[Metre, -1] :: HNil. Invert takes one input type and returns one output type:

trait Invert[A] {
	type Out
}
object Invert {
	type Aux[A, Out0] = Invert[A] { type Out = Out0 }
}

To inductively derive typeclass instances for inverting, we need to prove that:

1. We can derive an instance for HNil (the base case)
2. We can derive an instance for a non-empty HList (the n + 1 case), provided there is an existing instance for its tail (the n case)

The base case operates on HNil

implicit def baseCase: Invert.Aux[HNil, HNil] = new Invert[HNil] { type Out = HNil }

The inductive case assumes that the tail has an instance, and derives an instance for the head by negating the exponent:

implicit def inductiveCase[D, Exp <: XInt, NExp <: XInt, Tail <: HList,
  OutTail <: HList](
    implicit negateEv: OpInt.Aux[Negate[Exp], NExp],
    tailEv: Invert.Aux[Tail, OutTail]
): Invert.Aux[FieldType[D, Exp] :: Tail, FieldType[D, NExp] :: OutTail] =
  new Invert[FieldType[D, Exp] :: Tail] {
    type Out = FieldType[D, NExp] :: OutTail
}

When the compiler looks for the implicit instance for FieldType[Metre, 1] :: HNil:

• It finds that the inductiveCase method has a return type which fits the signature
• It can find the required evidence negateEv for negating 1 from singleton ops
• It requires evidence of an implicit instance for the tail HNil
• It finds that baseCase provides this evidence

So in hunting for implicit typeclass instance for the whole list, the compiler goes and finds instances for the tail (the n case), right up until the base. If we provide an inductive proof with a baseCase and an inductiveCase, we fit the bill for what the compiler needs.

Multiply

Now that we've tested a basic example of induction, we can go on to a more complex one.

We want to multiply two HLists of dimensions together. This means that the exponents should be added.

trait Multiply[A, B] {
  type Out
}
object Multiply {
  type Aux[A, B, Out0] = Multiply[A, B] { type Out = Out0 }
}

This is harder to make inductive because there are two input lists involved. Luckily, we only need to recurse over one of them, as we can pick dimensions from the other using shapeless's Selector. We will recurse over the left list and can pick elements from the right list.

Our base case can be the same:

implicit def baseCase: Multiply.Aux[HNil, HNil, HNil] = new Multiply[HNil, HNil] {
  type Out = HNil
}

We can define the inductive case using the following logic: 1. Pick the exponent in the right list corresponding to the head dimension in the left list 2. Add the left and right exponents together 3. Filter the term from the right list to get the remaining elements 4. Look for a typeclass instance for the left list tail and the remaining elements in the right list

implicit def inductiveCase[D, R <: HList, LExp <: XInt , RExp <: XInt,
  OutExp <: XInt, RTail <: HList, LTail <: HList, OutTail <: HList](
  implicit pickEv: Selector.Aux[R, D, RExp],
  addEv: OpInt.Aux[LExp + RExp, OutExp],
  filterEv: FilterNot.Aux[R, FieldType[D, RExp], RTail],
  tailEv: Multiply.Aux[LTail, RTail, OutTail]
): Multiply.Aux[FieldType[D, LExp] :: LTail, R, FieldType[D, OutExp] :: OutTail] =
  new Multiply[FieldType[D, LExp] :: LTail, R] {
    type Out = FieldType[D, OutExp] :: OutTail
}

When the compiler looks for an implicit instance of multiply for FieldType[Metre, -1] :: HNil and FieldType[Metre, 3] :: HNil:

• It finds that the inductiveCase has a return type which fits the signature
• Given that the head of the left list is Metre, it selects the exponent for Metre from the right list
• It can find the evidence addEv to add the exponents -1 and 3
• It filters Metre from the right list to get HNil
• It requires evidence of an instance for HNil and HNil
• This is provided by the base case

The compiler can now find instances of Multiply, as long as a dimension appears in both the left and right lists. This can be extended to when a dimension doesn't appear by writing a few more inductive cases.

Divide

The reason we went to the effort of writing Invert and Multiply was to divide. Dividing a numerator by a denominator is as simple as inverting the denominator and multiplying it by the numerator. We can write this in a single non-inductive instance:

implicit def divide[L <: HList, R <: HList, RInverted <: HList,
  Divided <: HList](
  implicit invertEv: Invert.Aux[R, RInverted],
  multiplyEv: Multiply.Aux[L, RInverted, Divided]
): Divide.Aux[L, R, Divided] = new Divide[L, R] {
  type Out = Divided
}

That's far simpler than the work we've done before - we're just building on the typeclasses we wrote to do this.

Automatically derived instances

We can now have compile time dimensional analysis without writing out divide instances for every combination of dimensions:

val rocketMass = Quantity[FieldType[Kilogram, 1] :: HNil](40.0)
val fuelMass = Quantity[FieldType[Kilogram, 1] :: HNil](10000.0)
val exhaust = Quantity[FieldType[Metre, 1] :: FieldType[Second, -1] :: HNil](1000000.0)
val distance = Quantity[FieldType[LightYear, 1] :: HNil](77.0)

val rocketSpeed = times(exhaust, log(divide(add(rocketMass, fuelMass), rocketMass).value))
val time: Quantity[FieldType[Year, 1] :: HNil] = divide(distance, rocketSpeed)
// error: type mismatch; found: FieldType[LightYear, 1] :: FieldType[Metre, -1] :: FieldType[Second, 1] :: HNil; required: FieldType[Year, 1] :: HNil

Yay! We've solved the problem! It looks a lot more verbose than what we started with, but we can tidy this up by using extension methods:

implicit final class DoubleOps(val d: Double) {
  def ly: Quantity[FieldType[LightYear,1] :: HNil] = Quantity(d)
  def kg: Quantity[FieldType[Kilogram, 1] :: HNil] = Quantity(d)
  def yr: Quantity[FieldType[Year, 1] :: HNil]     = Quantity(d)
  def mps: Quantity[FieldType[Metre, 1] :: FieldType[Second, -1] :: HNil] = Quantity(d)
  def c: Quantity[FieldType[LightYear, 1] :: FieldType[Year, -1] :: HNil] = Quantity(d)
}

We can also add symbolic infix operators for add, times and divide:

case class Quantity[A](value: Double) {
  def +(q: Quantity[A]): Quantity[A] = Quantity(value + q.value)
  def *(v: Double): Quantity[A] = Quantity(value + v)
  def /[B](q: Quantity[B])(implicit d: Divide[A, B]): Quantity[d.Out] = Quantity(value / q.value)
}

We've reached Libra!

We started with doubles:

val rocketSpeed = 1000000.0 * log((40.0 + 10000.0) / 40.0)
val time = 77.0 / rocketSpeed

And we finished with compile time dimensional analysis:

val rocketSpeed = 1000000.0.mps * log(((40.0.kg +10000.0.kg) /40.0.kg).value)
val speedConversion = 300000000.0.mps / 1.c
val speedInC = rocketSpeed / speedConversion
val time = 77.0.ly / speedInC
//time: Quantity[FieldType[Year, 1] :: HNil] = Quantity(4180.65274634)

The code isn't more verbose - if anything, it's more explanatory and just as easy to work with.

Rolling this out to more problems

All we need to provide for the business logic of our rocket launch problem are the dimensions and DoubleOps. We could roll this out to any other problem. Let's say we wanted to do a currency conversion between GBP and DKK:

val exchangeRate: Quantity[FieldType[DKK, 1] :: FieldType[GBP, -1] :: HNil] =
   currentExchangeRate()
Val krone: Quantity[FieldType[DKK, 1] :: HNil] = 10.gbp * exchangeRate

We get dimensional analysis for any problem domain out of the box!

Most of the code we've written is library code. In fact, it's Libra code! Libra is a dimensional analysis library based on typelevel induction. It performs compile time dimensional analysis to any problem domain. It also uses spire for its numeric typeclasses, so can be used for far more than just doubles.

Conclusion

It's been a long way from the humble Double. We started with basic types, explored dependent types, took a look at Typelevel Scala along the way, before finally ending up performing typelevel induction. As a result, we've managed to achieve compile time dimensional analysis for any problem. If you're curious about typelevel induction take a look at the Libra codebase for more examples. Enjoy!