Software Engineering -- An Idea Whose Time Has Gone?

That's a powerful statement over there. One I'd never dare state (aloud, at least :). They are not my words, though, they are Tony DeMarco's, who famously said "You can't control what you can't measure" in his book, Controlling Software Projects: Management, Measurement, and Estimation, Prentice Hall/Yourdon Press, 1982.

Well, he hasn't changed his mind about that, but he seems to have come to the conclusion that the relative value of that, particularly balanced against the costs of measuring -- which he once discounted -- is not what he once thought. You can read it on this IEEE-published article.

Getting Around Type Erasure with Manifests

A friend posed an interesting problem to me yesterday. He wanted a registry-like structure in which he would put pairs of key/value, and then get them back. However, he wished for it to be type-safe. Specifically, if he put a List[Int] inside it, and tried to get a List[String] out, he shouldn’t get anything.

This is generally a problem because of type erasure, a restriction of the JVM. At runtime, a List object knows it is a List, but not of what kind of element.

Alas, an obscure, and experimental, feature of Scala let you get around that. It’s the Manifests. A Manifest is class whose instances are objects representing types. Since these instances are objects, you can pass them around, store them, and generally call methods on them. With the support of implicit parameters, it becomes a very powerful tool:

object Registry {
  import scala.reflect.Manifest
  
  private var _map= Map.empty[Any,(Manifest[_], Any)] 
  
  def register[T](name: Any, item: T)(implicit m: Manifest[T]) {
    _map = _map(name) = (m, item)
  }
  
  def get[T](key:Any)(implicit m : Manifest[T]): Option[T] = {
    val o = _map.get(key)
   
    o match {
      case Some((om: Manifest[_], s : Any)) =>
        if(om <:< m) Some(s.asInstanceOf[T]) else None
      case _ => None
    }
  }
}

scala> Registry.register('a, List(1,2,3))

scala> Registry.get[List[Int]]('a)
res6: Option[List[Int]] = Some(List(1, 2, 3))

scala> Registry.get[List[String]]('a)
res7: Option[List[String]] = None

Above, Manifest is being passed as an implicit parameter – defined inside scala.reflect.Manifest – based on the type parameter T, and stored along the value. When you put a value in the register, it is not necessary to specify its type, as this gets inferred.

When you take a registry out, though, you specify the type the type you want (inference won’t help you here), and, if the key does not exist or does not store the desired type, it will return None. If a key exists and the type is correct, it will return Some(value).

Note that the operator <:< tests for subclassing. That is, the desired type m must be a superclass of the stored type om.

Delimited Continuations Explained (in Scala)

One of the promised features for Scala 2.8 is a delimited continuations plugin. A brief description was given in the Scala site, though the examples were, in my opinion, badly chosen.

Though my understanding is still pretty bare, I intend to shed a bit of light on it here. This is not directed at the experts, the people who will actually use them regularly, but at the people who will see them occasionally and will wonder what the heck is happening. Also, I might quite possibly have made mistakes below. If you do find mistakes, please correct me and I’ll fix the post.

First, let’s discuss why Scala is adopting continuations, and what, exactly, are them.

Continuations are, in essence, a program flow control structure, like conditional statements, loops, exceptions, etc. As mentioned in Scala’s Continuations Paper, classical, or full, continuations can be seen as a functional version of GOTO statements. If you haven’t heard of GOTO statements before, it’s enough to know they are essentially extinct for good reason.

Still quoting, delimited continuations are more like regular functions and less like GOTOs, and this makes them more powerful than the classical ones. So, basically, we are looking into a way of changing the execution flow of a program.

Using continuations is not easy. In fact, it’s difficult and error-prone, according to some references. But they can be very efficient performance-wise, and they allow one to build libraries with some interesting flow control abstractions built on top of continuations.

For these reasons, I doubt most Scala users will have much contact with them, and I suspect the ones who do will regret the necessity. The people writing libraries and frameworks, though, will likely be grateful that they are available. In that respect, continuation is not unlike some other advanced Scala features.

Ok, so let’s get to them. In Scala, continuations will be used through two keywords (or functions?), “reset” and “shift”. Like “catch” and “throw”, a “shift” always happen inside a “reset”, the later being responsible for delimiting the scope of the continuation. Before we proceed, let’s give some usage examples, which I’ll then try to explain.

Example 1:

reset {
shift { k: (Int=>Int) =>
 k(7)
} + 1
}
// result: 8

Example 2:

def foo() = {
println("Once here.")
shift((k: Int => Int) => k(k(k(7))))
}
def bar() = {
1+ foo()
}
def baz() = {
reset(bar() * 2)
}
baz()  // result 70

Example 3:

reset {
 println(1)
}
//prints: 1

reset {
 println(1)
 shift { (cont: Unit => Unit) => }
 println(2)
}
//prints: 1

reset {
 println(1)
 shift { (cont: Unit => Unit) =>
     println(2)
 }
 println(3)
}
//prints: 1 2

reset {
 println(1)
 shift { (cont: Unit => Unit) =>
     cont()
     println(2)
 }
 println(3)
}
//prints: 1 3 2

reset {
 println(1)
 shift { (cont: Unit => Unit) =>
     cont()
     cont()
     println(2)
 }
 println(3)
}
//prints: 1 3 3 2

So we have two examples of continuations as expressions, and one as purely flow control. The first example seems pretty simple: the whole shift is replaced with “7”, and, therefore, reset returns “7+1”. The third example shows that, in fact, the code before shift is executed once, the code after shift as many times as the function inside shift is repeated, and the code inside the shift itself is executed once, going back and forth to the code after it.

So we look at the second example, and see that… the code before shift seems to get executed many times! Actually, the method receiving the result of shift is being executed many times, once for each time “shift” return, but the statement right before shift only gets executed once.

These examples were probably enough to… make things more confusing, right? So let’s understand what is happening. Scala’s continuations are based on a selective Continuation Passing Style (CPS) transformation. The Wikipedia page on CPS is actually pretty good, so, please, take a moment to peruse it, and understand the examples given at the beginning, and then come back here.

Ready? Well, one thing that immediately came to mind is that this looks like programming backwards! They call it “inside-out” in that page, which is more appropriate, I suppose, but one has to pay attention to first impressions. :-)

What is happening there is that the innermost functions of normal programming style are outmost, and vice versa. You’ll see that this is what is happening with shift and reset.

I’ll now explain a bit of what is happening in a continuation. Notice, first, that the code associated with the shifts are functions: they are all in the form { k => … }. So, the first thing to keep in mind is that shifts are functions taking a parameter.

Now, pay attention, again, to the very first example. It defines the type of k as being “Int => Int”. In the first example, then, shift is a function which takes as parameter another function. This is actually true of all shifts, even if the function being passed to shift takes no parameters and returns Unit.

So a shift is a function which receives another function as parameter. As the body of the shift gets executed, that function-parameter may get called or not. From the second example, it is pretty obvious that the code after the shift only gets executed if that function-parameter is called.

By now you probably have an uneasy intuition that that parameter is, in fact, the code after the shift. That’s essentially right, but let’s further into this.

The result of reset is the result of the code inside shift. So, in the first example, “8”, the result of reset, is also the result of the shift code, “k: (Int => Int) => k(7)”. It should be obvious, then, that the value passed to “k” is “x: Int => x + 1”, which is a function of Int => Int, as required. Now, take a look at the two lines below:

shift { k: (Int=>Int) => k(7) } + 1
x                               + 1

So, the function being passed as parameter is formed by taking the whole code from the shift on, and replacing shift with “x”. In effect, when k(7) is executed, this is like replacing the whole shift with 7, and then continuing until the end of reset. At that point, though, the execution resumes inside the shift. In this example there is nothing more in there, so the result is 8.

What is actually happening is this:

((k: Int => Int) => k(7))((x: Int) => x + 1)

As said before, we are passing a function to another. You can paste that code inside REPL and it will work as expected. You may also try this code here:

((k: Int => Int) => k(k(k(7))))((x: Int) => (1 + x) * 2)

Are you getting a feeling now for it? So, whenever you see reset/shift, first think of the whole reset expression, replacing the shift inside it with an unknown variable – we’ll keep calling it x for now:

Reset { before shift { k => ... k(something) ... } after }
Reset { before x after }

Now replace Reset with a function declaration of one parameter:

Reset { before x after }
x => { before x after }

What type is “x” supposed to be? You can find it inside the shift, by looking into what is being passed to “k” – the “something” in the example. We can now rewrite the whole reset statement, placing first the shift code, then the reset code as parameter, like we did before:

( k => ... k(something) ... ) {x => { before x after }}

And this is pretty much what happens when you use continuations.

As for the second example, though the “1 + “ looks to be to the “left” of shift, it actually gets executed after shift returns, so it is part of the function being passed to shift. You find the same situation when making tail recursive functions. If you tried “1 + recurse()”, it would not be tail recursive, because “+” would only be executed after the call to “recurse” returned.

Pretty easy, right? Let’s try something with two shifts:

type Monadic[+U, C[_]] = {
def flatMap[V](f: U => C[V]): C[V]
}

class Reflective[+A, C[_]](xs: Monadic[A,C]) {
def reflect[B](): A @cps[C[B], C[B]] = {
 shift { k:(A => C[B]) =>
   xs.flatMap(k)
 }
}
}

implicit def reflective[A](xs:Iterable[A]) =
new Reflective[A,Iterable](xs)

reset {
val left = List("x","y","z")
val right = List(4,5,6)

List((left.reflect[Any], right.reflect[Any]))
}
// result: List[(String, Int)] = List((x,4), (x,5), (x,6), (y,4), (y,5), (y,6), (z,4), (z,5), (z,6)

Don’t hate me, it is in good cause. :-) If you can understand this… well, you can understand this. But I’m sure you’ll start to really get a grip in it.

First, let’s discard the irrelevant parts of the code above. Everything but the “reflect” definition and the reset code is irrelevant. They are there just to make “reflect” available to a List, it’s plain-old Enrich My Library pattern, though the use of a structural type is cool. :-)

Now, let’s replace left.reflect and right.reflect with the actual code getting executed – I’ll change k into k1 and k2 so as not to confuse both, so we are left only with the reset code:

reset {
val left = List("x","y","z")
val right = List(4,5,6)

List((shift { k1:(String => List[Any]) => left.flatMap(k1)}, shift { k2:(Int => List[Any]) => right.flatMap(k2)}))
}

Since there are two shifts, should we have two unknowns being passed to shift? If you look at the type signatures for k1 and k2, the answer is no. Each shift gets a single parameter.

You may think that the “Any” type means one list is formed by a single element, and the other by a tuple, but that’s not true. You can change both types to “List[(String,Int)]”, and it would still work. No, the type signatures are right, the trick is in how they are composed together.

If you think it through simply in terms of flow control, the first shift will return a value (“x”), and then the second shift will return a value (4), producing a list. The flow then returns to the last shift, which will return another value (5), and then another still (6). At this point, the flow returns once again to the first shift, which will produce it’s next value (“y”). So the trick is in how we compose the functions.

Let’s see, first, both shifts:

((k1: String => List[(String,Int)]) => left.flatMap(k1))
((k2: Int => List[(String,Int)]) => right.flatMap(k2))

These can’t and won’t be changed. The trick is how we compose them:

((k1: String => List[(String,Int)]) => left.flatMap(k1))
((x: String) => ((k2: Int => List[(String,Int)]) => right.flatMap(k2))
             ((y: Int) => List((x,y)))
)

So what happens is that the function which takes a String and returns a List of (String,Int) that gets passed to the first shift is composed of both the reset function and the second shift function.

The function passed to the second shift does not need an additional String parameter because the string is already bound.

Now you have seen both the control flow of reset/shift, and how to translate them into common functions to understand what values they will produce. I hope it will be of help when the time comes in which you need to understand these structures.

Note: many things have changed since I wrote this blog. Continuations have gone through minor changes in the annotations, and collections have gone through big changes. I haven't checked out all examples yet, but I know the one with the two lists doesn't work anymore. The following code will work for that example, though I wish I could make it more general than it presently is:

type Monadic[+U, C[_]] = {
  def flatMap[V, That](f: U => C[V])(implicit cbf: collection.generic.CanBuildFrom[C[U], V, That]): That
}

class Reflective[+A](xs: Monadic[A,Traversable]) {
  def reflect[B](): A @cps[Traversable[B]] = {
    shift { k:(A => Traversable[B]) =>
      xs.flatMap(k)(collection.breakOut) : Traversable[B]
    }
  }
}

implicit def reflective[A](xs:Traversable[A]) = new Reflective[A](xs)

reset {
  val left = List("x","y","z")
  val right = List(4,5,6)

  List((left.reflect[Any], right.reflect[Any]))
}

Matrices 6

In this last installment of the Matrices series, we’ll look into making the Matrix class parameterized. What does that means? Basically, we want to have Matrix declared like this:

class Matrix[T]

So that we can create matrices of Ints, Doubles, or whatever other numeric type we need. This is not a trivial exercise, because of the way numeric types are defined in Java, and Scala’s backward compatibility with it. Basically, there is no common superclass to them.

Paul Chisuano offered two alternatives around that. One would involve the “Enrich My Library” pattern, but I do not think it would offer any gains over MatrixInt as it was built, and I see a problem or two implementing it.

We’ll discuss the other solution he proposes, which is, indeed, very clever. If we were building this in Java, would idea would be creating an interface with the operations we need, and a concrete class for each type we need. For instance, let’s call our interface Scalar, and our concrete classes implementing it IntScalar, DoubleScalar, etc. We would then define a Matrix class parameterized with an “upper bound” of Scalar, created from an Array of Arrays of Scalar. That would work, but can be quite slow.

Instead of that, we’ll use one of the more wizardry of Scala’s features, the implicit type conversion argument. You have probably seen things like this before:

implicit def xToY(x: X): Y = Y(x)

This definition makes it possible to convert an object of class X to an object of class Y without explicitly doing so. If I ever need an object of class Y and I have an object of class X instead, Scala will apply that definition to object X without any explicit statement on my part to do so.

Even if you haven’t seen that before, you almost certainly have used it. See the two examples below:

scala> val a = 1.0; val b = 1
a: Double = 1.0
b: Int = 1

scala> a + b
res12: Double = 2.0

scala> val c = Map('a -> 1)
c: scala.collection.immutable.Map[Symbol,Int] = Map('a -> 1)

In each case, there is an implicit conversion at work. These definitions can, in fact, be found inside the object scala.Predef, which is imported automatically into scope of every Scala program. The first one is obvious, the second one, not so much. To understand the second, you have to realize that class Symbol does not define the method “->”, and that “->” is not a Scala keyword. I’ll leave as an exercise to find what implicit conversion is being used here.

What that means for us is that we can create a Scalar class, as well as implicit conversions from our desired types to it. Scalar would look like this:

abstract class Scalar[T](self: T) {
 def +(other: T): T
 def -(other: T): T
 def /(other: T): T
 def %(other: T): T
 def *(other: T): T
}

One thing to note here is that none of the operations defined for Scalar take another Scalar as argument, nor do they return a Scalar. The reason for this is that we do not want to store Scalars. They are bulky and inefficient. We just want to quickly convert an Int or Double into Scalar, and then produce, through an Int or Double operation, and Int or Double result to be stored in an Array, or returned as result of determinant or similar methods.

Even without implicits, we could define Matrix like this:

class Matrix[T](self: Array[Array[T]], toScalar: T => Scalar[T])

We could then use toScalar everywhere we do an operation over a type T. With implicit, we could just add an implicit definition to the top of class Matrix, so we don’t need bother explicitly using toScalar. But we don’t need to do even that. Look at the following definition:

class Matrix[T](self: Array[Array[T]])(implicit toScalar: T => Scalar[T])

This definition of Matrix constructor, with two sets of () for arguments, is said to be curried. I won’t be explaining what that means here, because what interests us is the implicit keyword in there. The important thing here is that when you define something like that, then two things useful for us happen:

1. If there is an implicit definition of the type “T => Scalar[T]” in scope when we call Matrix, then it will be automatically passed to the constructor. Or, in other words, you can call “new Matrix(Array(Array(1)))” anywhere there is an implicit definition from Int to Scalar[Int], without having to pass that implicit definition as well.
2. There is no need to define an implicit conversion inside Matrix itself. Any implicit conversion passed as an argument will already be used inside the class.

This rule is also valid for implicit conversions arguments on method definitions, by the way.

Anyway, this is so useful and powerful that Scala has another way of saying the same thing:

class Matrix[T <% Scalar[T]](self: Array[Array[T]])

The “<%” keyword means there should be an implicit conversion between the first type and the second in scope where Matrix constructor is called, or that one such conversion must be explicitly passed.

So, we search&replace every instance of MatrixInt for Matrix[T], define our constructor as above (plus minor details such as “private” and “val”), replace Int with T where appropriate (remember that row and column coordinates are still Ints!), create the Scalar class, the subclasses for the types we want, the implicit conversions, and we are done, right?

Not so. Scalar, alas, is not enough. If you look at the definitions of sign, isDivisibleBy and invert, you will see we mix “T” with -1, 0 and 1. Since we do not know that -1, 0 and 1 can be mixed with “T”, we still need something else. Alas, what we need is actually very simple: an implicit definition between Int and T, so that when we use -1, it can be converted to the proper type, be it Int, Double or whatever.

We have to be careful of one thing, though. See the two lines below:

map Function.tupled((value, col) => value * cofactor(0, col).determinant * sign(0, col))
map Function.tupled((value, col) => sign(0, col) * value * cofactor(0, col).determinant)

The first line works, because it is converted into this:

map Function.tupled((value, col) => toScalar(value) * toScalar(cofactor(0, col).determinant) * intToT(sign(0, col)))

The second line doesn’t work, as it would need to be converted into this:

map Function.tupled((value, col) => toScalar(intToT(sign(0, col))) * toScalar(value) * cofactor(0, col).determinant)

On the first line there was at most one implicit conversion for each parameter. On the second line, though, since T doesn’t have the method “*”, and Int doesn’t have the method “*(T)”, two implicit conversions would need to be applied.

Scala won’t apply two implicit conversions on the same element. In this case, we can make “sign” return T, so the extra implicit conversion will be applied inside “sign”, thus avoiding this rule. We need to be careful, however, about such situations. As it happens, MatrixInt, as it was written, won’t present any such problems.

Anyway, we need, thus, two implicit conversions. The “<%” syntax isn’t enough for our particular case, so we have to do it the other way. Here’s our Matrix declaration:

class Matrix[T] private (private val self: Array[Array[T]])
                       (implicit toScalar: T => Scalar[T], intToT: Int => T) {

Note that I write “implicit” only once in the declaration. All parameters in these last parentheses will be defined as implicit. Also, please note that implicit parameters must come last in such a declaration.

This does take care of almost everything, with a bit of judicious editing on various declarations to make the signatures match. There are some places, though, where that is not the case. Look at these lines inside the method “toString”:

   val maxNumSize = self.projection flatMap (_ map (_.toString.size)) reduceLeft (_ max _)
   val numFormat = "%"+maxNumSize+"d"
   def formatNumber(n: T) = numFormat format n

We obviously can’t use “%d” for Doubles, for one thing. Also, maxNumSize doesn’t quite cut it, since we would like numbers to be aligned on the decimal point (or comma, depending on your locale), so we need to know at least two sizes. And then there are complex numbers, and who knows what else, so we don’t really know what information might be needed.

Therefore, I feel it is better to leave the choice inside Scalar itself. I propose to add the method below to Scalar, define it for Int as shown next, and then use as show last:

abstract class Scalar[T](self: T) {
 ...
 def numFormat(xs: Iterator[T]): String
}

class IntScalar(self: Int) extends Scalar(self) {
 ...
 def numFormat(xs: Iterator[Int]) = {
   val maxNumSize = xs map (_.toString.size) reduceLeft (_ max _)
   "%"+maxNumSize+"d"
 }
}

 lazy val numFormat = this(0,0).numFormat(self flatMap (x => x) elements)
 override def toString = {
   def formatNumber(n: T) = numFormat format n
   ...

It would have been more appropriate to define that method inside a companion object, but this way is simpler. Note, as well, that “numFormat” became a visible property, so we can re-use all this logic inside LinearEquations.

Another place that needs fixing is inverse. We use the following rule:

 lazy val inverse: Option[Matrix[T]] =
   if (determinant == 0 || !isDivisibleBy(determinant))

Now, isDivisibleBy is implemented in terms of “%”, but that is just plain wrong for Double numbers. We could “fix” it by defining “%” to return 0, but that would make “%” act in a very unexpected way. So, instead, I’ll just define an “isDivisibleBy” operator inside Scalar:

 def isDivisibleBy(other: T): Boolean

The last problem, as far as I can see, is with the equality methods. You’ll get erasure warnings in them if you use “Matrix[T]” on the case statement of equals or inside asInstanceOf on canEqual. Instead, use “Matrix[_]”, to indicate you don’t care what type T is. Since we are doing a deepEquals, T will sort itself out.

The Scalar class and subclasses for Int and Double look like this:

abstract class Scalar[T](self: T) {
 def absolute: T
 def +(other: T): T
 def -(other: T): T
 def /(other: T): T
 def %(other: T): T
 def *(other: T): T
 def numFormat(xs: Iterator[T]): String
 override def numberSign: String
 def isDivisibleBy(other: T): Boolean
}

class IntScalar(self: Int) extends Scalar(self) {
 override def absolute = self.abs
 override def +(other: Int) = self + other
 override def -(other: Int) = self - other
 override def /(other: Int) = self / other
 override def %(other: Int) = self % other
 override def *(other: Int) = self * other
 override def numFormat(xs: Iterator[Int]) = {
   val maxNumSize = xs.map(_.toString.size).foldLeft(0)(_ max _)
   "%"+maxNumSize+"d"
 }
 override def numberSign = if (self < 0) "-" else "+"
 override def isDivisibleBy(other: Int) = (self % other) == 0
}

class DoubleScalar(self: Double) extends Scalar(self) {
 override def absolute = self.abs
 override def +(other: Double) = self + other
 override def -(other: Double) = self - other
 override def /(other: Double) = self / other
 override def %(other: Double) = self % other
 override def *(other: Double) = self * other
 override def numFormat(xs: Iterator[Double]) = {
   val numbers = xs.toList
   val maxIntSize = (
     numbers
     .map(_.toString.takeWhile(c => c != '.' && c != ',').length)
     .foldLeft(0)(_ max _)
   )
   val maxDecimalSize = (
     numbers
     .map(_.toString.dropWhile(c => c != '.' && c != ',').length)
     .foldLeft(0)(_ max _)
   ) - 1
   "%"+maxIntSize+"."+maxDecimalSize+"f"
 }
 override def numberSign = if (self < 0) "-" else "+"
 override def isDivisibleBy(other: Double) = true
}

Note that we transform “xs” into a List inside DoubleScalar. That’s because we are iterating through it twice, which isn’t possible with an Iterator. Also note we added “absolute” and “numberSign”, which are used by LinearEquations.

The method “absolute” is not called “abs”, though, because “abs” is defined on RichInt, RichDouble, etc. So if we defined it for IntScalar, and then tried to use “abs” on an Int, the compiler would not know whether to use RichInt’s abs or IntScala’s abs.

We’ll also add the following two lines to the Matrix companion-object:

 implicit def intToScalar(n: Int): Scalar[Int] = new IntScalar(n)
 implicit def doubleToScalar(n: Double): Scalar[Double] = new DoubleScalar(n)

This way these two types, Int and Double, will automatically become available after importing “Matrix._”. The user can arbitrarily extend Matrix to cover any type, though, by creating a subclass of Scalar and defining an implicit conversion from the desired type to it.

Performance-wise, the boxing of the numeric types for each operation is detrimental, but it can probably be optimized by either the compiler or the JIT. In Scala 2.8 there will also be the option of marking a class “specialized”, which is creates specially optimized versions of a class for Java’s “primitive” types. I’m not sure it applies in this situation, though.

I’ll finish with the revised code for Matrix. If you have any questions or suggestions, please send them to me and I’ll try to address them. I hope you found this series useful!

class Matrix[T] private (private val self: Array[Array[T]])
                       (implicit toScalar: T => Scalar[T], intToT: Int => T) {
 import Matrix._

 val rows = self.size
 val cols = self(0).size

 private def _row(n: Int): Array[T] = self(n)
 private def _col(n: Int): Array[T] = self map (_(n))

 def row(n: Int): Array[T] = _row(n) map (x => x)
 def col(n: Int): Array[T] = _col(n)

 def apply(i: Int, j: Int): T =
   if (i >= rows || j >= cols || i < 0 || j < 0)
     throw new IndexOutOfBoundsException
   else
     self(i)(j)

 def update(row: Int, col: Int, newValue: T): Matrix[T] =
   new Matrix(createArrays(rows, cols, (i,j) =>
     if (row == i && col == j) newValue else this(i,j)))

 def update(what: Dimension, where: Int, newValue: Array[T]): Matrix[T] =
   what match {
     case Row => replaceRow(where, newValue)
     case Column => replaceCol(where, newValue)
   }
    
 def replaceCol(col: Int, newValue: Array[T]): Matrix[T] =
   new Matrix(createArrays(rows, cols, (i,j) =>
     if (col == j) newValue(i) else this(i,j)))

 def replaceRow(row: Int, newValue: Array[T]): Matrix[T] =
   new Matrix(createArrays(rows, cols, (i,j) =>
     if (row == i) newValue(j) else this(i,j)))

 override def equals(other: Any): Boolean = other match {
   case that: Matrix[_] =>
     that.canEqual(this) && self.deepEquals(that.self)
   case _ => false
 }

 def canEqual(other: Any): Boolean = other.isInstanceOf[Matrix[_]]
  
 def *(n: T): Matrix[T] =
   new Matrix(createArrays(rows, cols, (row, col) => this(row,col) * n))

 def /(n: T): Matrix[T] =
   new Matrix(createArrays(rows, cols, (row, col) => this(row,col) / n))

 def %(n: T): Matrix[T] =
   new Matrix(createArrays(rows, cols, (row, col) => this(row,col) % n))

 def +(other: Matrix[T]): Matrix[T] =
   if (rows != other.rows || cols != other.cols)
     throw new IllegalArgumentException
   else
     new Matrix(createArrays(rows, cols, (row, col) => this(row,col) + other(row,col)))

 def -(other: Matrix[T]): Matrix[T] =
   if (rows != other.rows || cols != other.cols)
     throw new IllegalArgumentException
   else
     new Matrix(createArrays(rows, cols, (row, col) => this(row,col) - other(row, col)))

 def *(other: Matrix[T]): Matrix[T] =
   if (cols != other.rows)
     throw new IllegalArgumentException
   else
     new Matrix(createArrays(rows, other.cols, (row, col) =>
       this._row(row) zip other._col(col) map Function.tupled(_*_) reduceLeft (_+_)
     ))

 def **(n: Int): Matrix[T] =
   if (rows != cols)
     throw new UnsupportedOperationException
   else n match {
     case 0 => unitMatrix(rows)
     case 1 => this
     case 2 => this * this
     case negative if negative < 0 =>
       (this ** negative.abs).inverse getOrElse (throw new UnsupportedOperationException)
     case odd if odd % 2 == 1 => this ** (odd - 1) * this
     case even => this ** (even / 2) ** 2
   }

 def toArray = Matrix.clone(self)

 def transpose = new Matrix(createArrays(cols, rows, (row,col) => this(col,row)))

 def cofactor(row: Int, col: Int) = new Matrix(createArrays(rows-1, cols-1, (i,j) =>
   this(i + (if (i < row) 0 else 1), j + (if (j < col) 0 else 1))))

 protected def sign(row: Int, col: Int): T = if ((col + row) % 2 == 0) 1 else -1

 lazy val determinant: T =
   if (rows != cols)
     throw new UnsupportedOperationException
   else rows match {
     case 1 => this(0,0)
     case 2 => this(0,0)*this(1,1) - this(1,0)*this(0,1)
     case n => (
       _row(0).zipWithIndex
       map Function.tupled((value, col) => value * cofactor(0, col).determinant * sign(0, col))
       reduceLeft (_+_)
       )
     }

 def minor(row: Int, col: Int) = cofactor(row, col).determinant

 lazy val adjugate: Matrix[T] = new Matrix(createArrays(cols, rows, (row, col) =>
   minor(col, row) * sign(col, row)
 ))

 def isDivisibleBy(n: T): Boolean =
   self flatMap (row => row) forall (_ isDivisibleBy n)

 lazy val inverse: Option[Matrix[T]] =
   if (determinant == 0 || !isDivisibleBy(determinant))
     None
   else
     Some(adjugate / determinant)

 lazy val numFormat = this(0,0).numFormat(self flatMap (x => x) elements)

 override def toString = {
   val numFormat = this(0,0).numFormat(self flatMap (x => x) elements)
   def formatNumber(n: T) = numFormat format n
   val top = rows match {
     case 1 => _row(0) map formatNumber mkString ("< ", " ", " >")
     case _ => _row(0) map formatNumber mkString ("/ ", " ", " \\") // Fix highlighter: "
   }
   val middle = rows match {
     case 1 | 2 => Nil
     case _ => self.toList.tail.init map (_ map formatNumber mkString("| "," "," |"))
   }
   val bottom = rows match {
     case 1 => Nil
     case _ => List(_row(rows - 1) map formatNumber mkString ("\\ ", " ", " /"))
   }
  
   top :: middle ::: bottom mkString "\n"
 }
}

object Matrix {
 import Matrix._
 implicit def toMatrix[T](m : Array[Array[T]])
                         (implicit toScalar: T => Scalar[T], intToT: Int => T) = Matrix(m)

 def apply[T](array: Array[Array[T]])
             (implicit toScalar: T => Scalar[T], intToT: Int => T): Matrix[T] =
   new Matrix(clone(array))
 def apply[T](rows: Int, cols: Int)
             (implicit toScalar: T => Scalar[T], intToT: Int => T): Matrix[T] =
   Matrix(rows, cols, 0)
 def apply[T](rows: Int, cols: Int, value: T)
             (implicit toScalar: T => Scalar[T], intToT: Int => T): Matrix[T] =
   new Matrix(createArrays(rows, cols, ((_,_) => value)))
 def apply[T](rows: Int, cols: Int, f: (Int,Int) => T)
             (implicit toScalar: T => Scalar[T], intToT: Int => T): Matrix[T] =
   new Matrix(createArrays(rows, cols, f))

 def unitMatrix[T](n: Int)(implicit toScalar: T => Scalar[T], intToT: Int => T) =
   Matrix[T](n, n, (row: Int,col: Int) => (if (row == col) intToT(1) else intToT(0)))

 protected def createArrays[T](rows: Int, cols: Int, f: (Int, Int) => T) =
   for((i: Int) <- (0 until rows).toArray)
     yield for((j: Int) <- (0 until cols).toArray)
       yield f(i,j)

 protected def clone[T](a: Array[Array[T]]) =
   createArrays(a.size, a(0).size, (row, col) => a(row)(col))
  
 sealed abstract class Dimension
 case object Row extends Dimension
 case object Column extends Dimension

 implicit def intToScalar(n: Int): Scalar[Int] = new IntScalar(n)
 implicit def doubleToScalar(n: Double): Scalar[Double] = new DoubleScalar(n)
}

class LinearEquations[T](val m: Matrix[T], val r: Matrix[T])
                       (implicit toScalar: T => Scalar[T], intToT: Int => T) {
 require(m.rows == r.rows && r.cols == 1)

 def maxColumnSize(matrix: Matrix[T], isAbs: Boolean): Int = (
   matrix.toArray
   flatMap (row => row)
   map (value => (if (isAbs) value.absolute else value).toString.length)
   reduceLeft (_ max _)
 )

 lazy val maxColsSize = m.cols.toString.length

 def numberSign(n: T) = " %s " format n.numberSign
 def formatConstant(n: T) = r.numFormat format n
 def formatCoefficient(n: T) = m.numFormat format n.absolute
 def formatUnknown(index: Int) = "x("+("%0"+maxColsSize+"d" format index)+")"

 // On Scala 2.8, replace drop(1) with tail and first with head below
 def formatLine(coefficients: Array[T],
                constant: T,
                printUnknown: Int => String) =
   coefficients
   .zipWithIndex
   .drop(1)
   .foldLeft(formatCoefficient(coefficients.first)+printUnknown(0))(
     (x,y) => x+numberSign(y._1)+formatCoefficient(y._1)+printUnknown(y._2)
   )+" = "+formatConstant(constant)

 def makeString(printUnknown: Int => String) = (
   for(row <- 0 until m.rows)
   yield formatLine(m.row(row), r(row,0), printUnknown)
 ).mkString("\n")
  
 override def toString = makeString(formatUnknown)

 def getUnknowns(u: Matrix[T]): Int => String =
   (index: Int) => " * "+(u.numFormat format u(index,0))

 lazy val solutionInverse: Option[Matrix[T]] =
   if (m.inverse == None)
     None
   else
     Some(m.inverse.get * r)

 def solveInverse =
   if (solutionInverse == None)
     "There is no unique solution"
   else
     makeString(getUnknowns(solutionInverse.get))

 lazy val solutionCramer: Option[Matrix[T]] = {
   val vector = r.toArray.flatMap(x => x)
   if (m.determinant == 0)
     None
   else
     Some(Matrix(
       for((col: Int) <- (0 until m.cols).toArray)
       yield for(row <- Array(0)) // Array[T](...) doesn't work, as it requires T <: AnyRef
         yield m.replaceCol(col, vector).determinant / m.determinant
     ))
 }
    
 def solveCramer =
   if (solutionCramer == None)
     "There is no unique solution"
   else
     makeString(getUnknowns(solutionCramer.get))
}

Matrices 5b

To finish my last post I wrote a quick&dirty class to show MatrixInt in use. Alas, it was so dirty I couldn't stand the shame, so I rewrote it. Have fun. :-)

class LinearEquations(val m: MatrixInt, val r: MatrixInt) {
  require(m.rows == r.rows && r.cols == 1)
  
  def maxColumnSize(matrix: MatrixInt, isAbs: Boolean): Int = (
    matrix.toArray 
    flatMap (row => row) 
    map (value => (if (isAbs) value.abs else value).toString.length)
    reduceLeft (_ max _)
  )
  
  lazy val maxCoefficientsColumnSize =  maxColumnSize(m, true)
  lazy val maxConstantsColumnSize = maxColumnSize(r, false)
  lazy val maxColsSize = m.cols.toString.length
  
  def numberSign(n: Int) = if (n < 0) " - " else " + "
  def formatNumber(n: Int, size: Int) = "%"+size+"d" format n
  def formatConstant(n: Int) = formatNumber(n, maxConstantsColumnSize)
  def formatCoefficient(n: Int) = formatNumber(n.abs, maxCoefficientsColumnSize)
  def formatUnknown(index: Int) = "x("+("%0"+maxColsSize+"d" format index)+")"
  
  // On Scala 2.8, replace drop(1) with tail and first with head below
  def formatLine(coefficients: Array[Int], 
                 constant: Int, 
                 printUnknown: Int => String) =
    coefficients
    .zipWithIndex
    .drop(1)
    .foldLeft(formatCoefficient(coefficients.first)+printUnknown(0))(
      (x,y) => x+numberSign(y._1)+formatCoefficient(y._1)+printUnknown(y._2)
    )+" = "+formatConstant(constant)

  def makeString(printUnknown: Int => String) = (
    for(row <- 0 until m.rows)
    yield formatLine(m.row(row), r(row,0), printUnknown)
  ).mkString("\n")
    
  override def toString = makeString(formatUnknown)
  
  def getUnknowns(u: MatrixInt): Int => String =
    (index: Int) => " * "+u(index,0)

  lazy val solutionInverse: Option[MatrixInt] =
    if (m.inverse == None)
      None
    else
      Some(m.inverse.get * r) 

  def solveInverse =
    if (solutionInverse == None)
      "There is no unique solution"
    else
      makeString(getUnknowns(solutionInverse.get)) 

  lazy val solutionCramer: Option[MatrixInt] = {
    val vector = r.toArray.flatMap(x => x)
    if (m.determinant == 0)
      None
    else
      Some(MatrixInt(
        for((col: Int) <- (0 until m.cols).toArray)
        yield Array(
            m.replaceCol(col, vector).determinant / m.determinant 
          )
      ))
  }
      
  def solveCramer =
    if (solutionCramer == None)
      "There is no unique solution"
    else
      makeString(getUnknowns(solutionCramer.get))
}