by Zainab Ali on Jun 13, 2017
technical
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.
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.
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.
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:
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
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.
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.
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.
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.
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.
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.
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.
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.
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:
This is similar to the mathematical method of proof by induction.
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:
HNil (the base 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:
inductiveCase method has a return type which fits the signaturenegateEv for negating 1 from singleton opsHNilbaseCase provides this evidenceSo 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.
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:
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:
inductiveCase has a return type which fits the signatureMetre, it selects the exponent for Metre from the right listaddEv to add the exponents -1 and 3Metre from the right list to get HNilHNil and HNilThe 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.
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.
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 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.
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.
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!
