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