Showing posts with label matrix. Show all posts
Showing posts with label matrix. Show all posts

Thursday, July 2, 2009

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

Wednesday, July 1, 2009

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

Tuesday, June 30, 2009

Matrices 5

Our Matix class is now almost complete. There are few methods I intend to add to it, and one I intend to change. Most of this is pretty straight-forward, though, so I’ll be quick about it.

First, matrix transposition:

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

Pretty trivial. Next, determinant. To compute the determinant of a matrix, though, I need to compute some co-factors of it. The co-factor method is the only one I’m unhappy about. My original implementation was something like this:

for((i: Int) <- (0 until rows).toArray; if i != row)
yield for((j: Int) <- (0 until cols).toArray; if j != col)
yield f(i,j)

That’s a pretty good implemention, but then I went and factored that code out of it and into createArrays. Unfortunately, there is no “if” in createArrays, and I’m not willing to add them for the sake one one method, as this can have some performance impact, and createArrays is used everywhere.

I could create a “createArraysCond” or something. That would do, but there is only one method that would make use of it. I’m wary of premature refactoring. :-)

So I ended with the method below. Also, the determinant:

def cofactor(row: Int, col: Int) = new MatrixInt(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) = if ((col + row) % 2 == 0) 1 else -1

lazy val determinant: Int =
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 (_+_)
)
}

I factored out “sign” because it is used in another method, and makes the code much clearer. Since our matrix is immutable, I made determinant “lazy val”. It’s a val because I need compute it only once – it’s value won’t change. It’s lazy because, that way, it won’t get computed until someone actually uses it. Otherwise, all matrices would have their determinant computed at creation time. This way, if you never need the determinant of a matrix, it won’t ever get computed.

The last methods I want to create for Matrix are “adjugate”, which will compute the adjugate, or adjoint, matrix, and “inverse”, which will compute the inverse of a Matrix. Since we are dealing with integers, computing the inverse will be a bit more difficult than if we were using doubles. Anyway, let’s do it. First, let’s define a “/” operator between matrix and integer. We’ll also define a “%” operator, which is pretty useful when dividing integers.

def /(n: Int): MatrixInt =
new MatrixInt(createArrays(rows, cols, (row, col) => this(row,col) / n))

def %(n: Int): MatrixInt =
new MatrixInt(createArrays(rows, cols, (row, col) => this(row,col) % n))

Next, I’ll do the adjugate:

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

lazy val adjugate: MatrixInt = new MatrixInt(createArrays(cols, rows, (row, col) =>
minor(col, row) * sign(col, row)
))

This is a bit confusing because I’m inverting at the same time I’m computing the minors. This saves a transposition step, though.

Finally, I’ll compute the inverse matrix, after asserting the matrix does have an integer inverse:

def isDivisibleBy(n: Int): Boolean =
(this % n).self flatMap (row => row) forall (_ == 0)

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

We finally have something new to look at. In the method isDivisibleBy, we first use the method “flapMap” to unnest the arrays, producing an Array[Int] containing all values of the matrix. Then we use “forall” in that. The method “forall” returns true if, and only if, the function passed to it (“_ == 0”) is true for all elements of the collection on which “forall” was invoked.

Or, in other words, “forall” will check that (this % n) results in a matrix of zeros.

The other interesting thing here is the use of Option. If you haven’t seen Option before, you probably haven’t seen much Scala code before. Option is Scala’s way of saying “this is not available”. An Option is a class, with a subclass called “Some”, and a singleton object called “None”. So where, in other languages, you might see a “null” value, in Scala the library will usually return a None instead.

On the other hand, instead of returning the value you want, it returns Some(value). This way, you need to test to see if the value you want is actually available, and then get it. There are many ways of doing that, and it pays to look up the API for Option.

The very last thing I’ll do with MatrixInt is change the exponentiation to be able to receive negative values, and we’ll see here one way of getting the value for an Option:

def **(n: Int): MatrixInt =
if (rows != cols)
throw new UnsupportedOperationException
else n match {
case 0 => unitMatrixInt(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
}

We’ll finish this installment with a couple of methods to solve systems of linear equations. It is tempting to create the methods inside MatrixInt itself, but solving linear equations is something you can do using matrices, not a matrix operation in itself.

I’m not sure what a good abstraction would be. We might just create the functions inside an object, and import them. Instead, I’ll create a small class to do that. First, let’s get a toString method. This method is quick and dirty, and could use some improvements. In particular, we could align the columns, like MatrixInt does, and we could treat negative numbers to avoid showing things like “+ -5”. That, I’ll leave as an exercise! :-)

class LinearEquations(val m: MatrixInt, val r: MatrixInt) {
require(m.rows == r.rows && r.cols == 1)

override def toString =
(for (rowNumber <- 0 until m.rows;
row = m.row(rowNumber).zipWithIndex)
yield (row
map Function.tupled(_+"x("+_+")")
mkString (""," + "," = ")) + r(rowNumber,0))
.mkString("\n")
}

To this class, we’ll add two solvers. First, through the inverted matrix, which should be pretty easy. We’ll also add a method which prints the equation with the unknowns replaced with the solution, so we can easily verify it by eye:

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
(for (rowNumber <- 0 until m.rows;
row = m.row(rowNumber).zipWithIndex)
yield (row
map Function.tupled((value, col) => value+"*"+solutionInverse.get(col,0))
mkString (""," + "," = ")) + r(rowNumber,0))
.mkString("\n")

We can already see the potential for refactoring. Now, the Cramer’s Rule. This is more difficult, because we need to compute each element of the solution. Though a bit more complex, it isn’t terribly so:

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
(for (rowNumber <- 0 until m.rows;
row = m.row(rowNumber).zipWithIndex)
yield (row
map Function.tupled((value, col) => value+"*"+solutionCramer.get(col,0))
mkString (""," + "," = ")) + r(rowNumber,0))
.mkString("\n")

More refactoring potential. We’ll finish this installment with an example of LinearEquations in action, followed by the complete source code for what we have done so far.

scala> m
res5: Array[Array[Int]] = Array(Array(5, -4), Array(6, -5))

scala> r
res6: Array[Array[Int]] = Array(Array(2), Array(1))

scala> new LinearEquations(m, r)
res2: LinearEquations =
5x(0) + -4x(1) = 2
6x(0) + -5x(1) = 1

scala> res2.solveInverse
res3: java.lang.String =
5*6 + -4*7 = 2
6*6 + -5*7 = 1

scala> res2.solveCramer
res4: java.lang.String =
5*6 + -4*7 = 2
6*6 + -5*7 = 1

And now, the source code:

package com.blogspot.dcsobral.matrix

class MatrixInt private (private val self: Array[Array[Int]]) {
import MatrixInt._

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

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

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

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

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

def update(what: Dimension, where: Int, newValue: Array[Int]): MatrixInt =
what match {
case Row => replaceRow(where, newValue)
case Column => replaceCol(where, newValue)
}

def replaceCol(col: Int, newValue: Array[Int]): MatrixInt =
new MatrixInt(createArrays(rows, cols, (i,j) =>
if (col == j) newValue(i) else this(i,j)))

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

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

def canEqual(other: Any): Boolean = other.isInstanceOf[MatrixInt]

def *(n: Int): MatrixInt =
new MatrixInt(createArrays(rows, cols, (row, col) => this(row,col) * n))

def /(n: Int): MatrixInt =
new MatrixInt(createArrays(rows, cols, (row, col) => this(row,col) / n))

def %(n: Int): MatrixInt =
new MatrixInt(createArrays(rows, cols, (row, col) => this(row,col) % n))

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

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

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

def **(n: Int): MatrixInt =
if (rows != cols)
throw new UnsupportedOperationException
else n match {
case 0 => unitMatrixInt(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 = MatrixInt.clone(self)

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

def cofactor(row: Int, col: Int) = new MatrixInt(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) = if ((col + row) % 2 == 0) 1 else -1

lazy val determinant: Int =
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: MatrixInt = new MatrixInt(createArrays(cols, rows, (row, col) =>
minor(col, row) * sign(col, row)
))

def isDivisibleBy(n: Int): Boolean =
(this % n).self flatMap (row => row) forall (_ == 0)

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

override def toString = {
val maxNumSize = self.projection flatMap (_ map (_.toString.size)) reduceLeft (_ max _)
val numFormat = "%"+maxNumSize+"d"
def formatNumber(n: Int) = 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 MatrixInt {
import MatrixInt._
implicit def toMatrixInt(m : Array[Array[Int]]) = MatrixInt(m)

def apply(array: Array[Array[Int]]): MatrixInt = new MatrixInt(clone(array))
def apply(rows: Int, cols: Int): MatrixInt = MatrixInt(rows, cols, 0)
def apply(rows: Int, cols: Int, value: Int): MatrixInt =
new MatrixInt(createArrays(rows, cols, ((_,_) => value)))
def apply(rows: Int, cols: Int, f: (Int,Int) => Int): MatrixInt =
new MatrixInt(createArrays(rows, cols, f))

def unitMatrixInt(n: Int) = MatrixInt(n, n, (row,col) => (if (row == col) 1 else 0))

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

protected def clone(a: Array[Array[Int]]) =
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
}

class LinearEquations(val m: MatrixInt, val r: MatrixInt) {
require(m.rows == r.rows && r.cols == 1)

override def toString =
(for (rowNumber <- 0 until m.rows;
row = m.row(rowNumber).zipWithIndex)
yield (row
map Function.tupled(_+"x("+_+")")
mkString (""," + "," = ")) + r(rowNumber,0))
.mkString("\n")

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
(for (rowNumber <- 0 until m.rows;
row = m.row(rowNumber).zipWithIndex)
yield (row
map Function.tupled((value, col) => value+"*"+solutionInverse.get(col,0))
mkString (""," + "," = ")) + r(rowNumber,0))
.mkString("\n")

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
(for (rowNumber <- 0 until m.rows;
row = m.row(rowNumber).zipWithIndex)
yield (row
map Function.tupled((value, col) => value+"*"+solutionCramer.get(col,0))
mkString (""," + "," = ")) + r(rowNumber,0))
.mkString("\n")
}

Sunday, June 28, 2009

Matrices 4

We have a reasonable Matrix class at this point, but it still has some deficiencies. There are a few common matrix operations we can add to it, and there are some problems related to the usability of the class. Let’s see them.

Equality

If you try “unitMatrixInt(2) == unitMatrixInt(2)”, you’ll get false as result. This happens because we did not provide an “equals” method, and the default “equals” method compares object references. Or, in other words, it tests if the objects are the same, not if the objects are equal. If that’s still unclear, see the following example:

scala> unitMatrixInt(2) == unitMatrixInt(2)
res8: Boolean = false

scala> val a = unitMatrixInt(2)
a: teste.MatrixInt =
/ 1 0 \
\ 0 1 /

scala> val b = unitMatrixInt(2)
b: teste.MatrixInt =
/ 1 0 \
\ 0 1 /

scala> val c = a
c: teste.MatrixInt =
/ 1 0 \
\ 0 1 /

scala> a == a
res9: Boolean = true

scala> a == b
res10: Boolean = false

scala> a == c
res11: Boolean = true

Writing an equals method can be tricky, but I’m not going to discuss here why such a method must be written this or that way, as there are good references on the web. I usually recommend this one.

I will talk about the syntactical aspect of making such a method. First, though the examples show the method “==”, I keep talking about “equals”. That’s because “==”, as well as “!=”, in Scala is implemented in terms of the “equals” method. So if you want to change them, you implement “equals”.

Now, it’s not enough just to call it “equals”. It has to have the correct signature, which is to receive a single parameter of type Any, and return a Boolean. The method below will be good enough, with the addition of a canEqual method, and a small change in our class declaration:

class MatrixInt private (private val self: Array[Array[Int]]) {
...
override def equals(other: Any): Boolean = other match {
case that: MatrixInt =>
that.canEqual(this) && self.deepEquals(that.self)
case _ => false
}

def canEqual(other: Any): Boolean = other.isInstanceOf[MatrixInt]

We changed the constructor so that self got declared as a “private val”, otherwise it won’t be visible even to other objects of the same class. We could also get a copy of a MatrixInt’s array representation, and compare against that, but that can make equals quite slow. To get such a copy, though, we would use the “toArray” method, which I thought I had presented in the first installment of these posts, but I now find I didn’t. Let’s correct that oversight:

def toArray = MatrixInt.clone(self)

We produce a copy of the array, so that changes to it won’t affect this object. Anyway, now the test works as expected:

scala> val a = unitMatrixInt(2)
a: teste2.MatrixInt =
/ 1 0 \
\ 0 1 /

scala> val b = unitMatrixInt(2)
b: teste2.MatrixInt =
/ 1 0 \
\ 0 1 /

scala> val c = a
c: teste2.MatrixInt =
/ 1 0 \
\ 0 1 /

scala> a == a
res12: Boolean = true

scala> a == b
res13: Boolean = true

scala> a == c
res14: Boolean = true


Updating

We also want to check to see it returns false for different matrices. To do that, though, we would like to produce a different matrix easily. Unfortunately, we have no way to do that. Well, we can generate unit matrices of different sizes, but if we just wanted to change a small value in one we have, we can do that.

Since our class is immutable, we do not expect it to actually change to something else, but when we do “1+1”, we do not change “1”, we just produce a new value. We can do the same thing, and enable individual cell changes in our matrix, like this:

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

Here are some tests, then:

scala> val a = unitMatrixInt(2)
a: teste3.MatrixInt =
/ 1 0 \
\ 0 1 /

scala> a == unitMatrixInt(3)
res20: Boolean = false

scala> a == unitMatrixInt(2)
res21: Boolean = true

scala> a == (a + a)
res22: Boolean = false

scala> a == a*2
res23: Boolean = false

scala> a == a * a // True! 1 == 1 * 1
res24: Boolean = true

scala> val b = a(1,0) = 1
b: teste3.MatrixInt =
/ 1 0 \
\ 1 1 /

scala> a
res25: teste3.MatrixInt =
/ 1 0 \
\ 0 1 /

scala> a == b
res26: Boolean = false

So, where’s “update”? It is actually used when declaring “b”. When you write “id(…) = value” in Scala, the compiler translates that to “id.update(…,value)”. So we could have also written this:

scala> val b = a.update(1,0,1)
b: teste3.MatrixInt =
/ 1 0 \
\ 1 1 /

Let’s also make it easy to change a whole column or a whole row. We can already do “a.row(0)” or “c.col(1)”, so doing “a.row(0) = …” would be nice. Scala has something similar, in that any method called “field_=” will be called if you write “id.field = …”. Unfortunately, if we try that here it won’t work, because “id.field” is already the name of a method which takes a parameter (the row or column). If it did not take any parameters, it would work. As it is, we’ll just have to make do with something else:

def replaceCol(col: Int, newValue: Array[Int]): MatrixInt =
new MatrixInt(createArrays(rows, cols, (i,j) =>
if (col == j) newValue(i) else this(i,j)))

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

And we can use it like this:

scala> val a = unitMatrixInt(3)
a: teste6.MatrixInt =
/ 1 0 0 \
| 0 1 0 |
\ 0 0 1 /

scala> a.replaceCol(0, a.col(2))
res39: teste6.MatrixInt =
/ 0 0 0 \
| 0 1 0 |
\ 1 0 1 /

There is one alternative, though, which I show below. I’m using the update method, and a symbol to identify what is being updated:

def update(what: Symbol, where: Int, newValue: Array[Int]): MatrixInt =
what match {
case 'row => replaceRow(where, newValue)
case 'col | 'column => replaceCol(where, newValue)
case _ => throw new IllegalArgumentException
}

scala> val a = unitMatrixInt(3)
a: teste7.MatrixInt =
/ 1 0 0 \
| 0 1 0 |
\ 0 0 1 /

scala> a('row, 1) = a.row(0)
res41: teste7.MatrixInt =
/ 1 0 0 \
| 1 0 0 |
\ 0 0 1 /

scala> a('column, 2) = a.row(0)
res42: teste7.MatrixInt =
/ 1 0 1 \
| 0 1 0 |
\ 0 0 0 /

Symbols are similar to String. The technical differences aren’t really much important – it’s their intent usage that differs. While strings are made to be shown to a user or received from them, a symbol is a piece of information exchanged between different parts of a program. They are very popular in some languages, but not so much on Scala.

One reason for their lack of popularity is that you don’t get any compile-time checks using them. If I misspelled “colunm” instead of “column”, the program would compile, and I’d get an exception at run-time – if I was lucky enough to catch it before the user.

Since we are using a statically typed language, we would like to catch such errors at compile-time. We can do that with Enumerations. Using a Scala enumeration is a bit more complex than other languages, as it’s not built in the language, but a library like anything else. It will provide just the checks we want, though.

First, we need to declare an object to be our enumeration. We’ll do this in the MatrixInt object-companion. We’ll also declare a type inside the object, to be used in our method declaration. It goes like this:

object Dimension extends Enumeration {
type Dimension = Value
val Row, Column = Value
}

“Value”, here, is both a subclass of Enumeration and a method which will return a new instance of type Value each time it is called. These instances will always be different between them. Note that when you declare “val Row, Column = Value”, Scala converts this into “val Row = Value; val Column = Value”, so “Value” gets invoked twice. If “Value” was called only once, then Row would be equal to Column, which is not something we want.

We might like to import the contents of the Dimension object inside the class too, so that we can write “what: Dimension” instead of “what: Dimension.Dimension” (or even “what: MatrixInt.Dimension.Dimension”, if we hadn’t imported object MatrixInt’s contents already). Simple line:

import MatrixInt.Dimension._

Now we can rewrite our method:

def update(what: Dimension, where: Int, newValue: Array[Int]): MatrixInt =
what match {
case Row => replaceRow(where, newValue)
case Column => replaceCol(where, newValue)
}

It’s important that Row and Column were declared capitalized – the first letter is uppercase. If we used “row” and “column”, then the match would simply assign “what” to “row”, like it did to “even” and “odd” in “**”, and execute the first statement. Since we used capitalized identifiers, the case statement compares “what” to the value of Row, as if it was a literal, and assigns nothing.

Note also that we dropped the default case. One has to be careful, though, because the compiler won’t warn against missing case statements. Here are some usage examples:

scala> a(Row, 2) = a.row(1)
<console>:13: error: not found: value Row
a(Row, 2) = a.row(1)
^

scala> a(Dimension.Row, 2) = a.row(1)
res1: teste.MatrixInt =
/ 1 0 0 \
| 0 1 0 |
\ 0 1 0 /

scala> import Dimension._
import Dimension._

scala> a(Row, 2) = a.row(1)
res2: teste.MatrixInt =
/ 1 0 0 \
| 0 1 0 |
\ 0 1 0 /

Can we do something about not getting warnings about missing case statements? As a matter of fact, yes, we can. Let’s rewrite once again our method. This time, we’ll use a case objects. Let’s show it first, and then discuss it. First, we replace the whole Dimension object by this:

sealed abstract class Dimension
case object Row extends Dimension
case object Column extends Dimension

Since we don’t have any nested declaration, we can do away with the extra import line. As for the rest… the update method doesn’t actually need any changes! The type, Dimension, has been declared. Row and Column are singleton objects, capitalized so that they can be used inside case statements as if they were literals.

There are two new things here. The “sealed” keyword, and the “case” keyword. Anyone reading about Scala has probably seen case classes. A case class is a class mostly like any other, but with a few methods automatically created by the compiler. A case object is similar to case class, and used instead of them when a case class would not declare any parameter in its constructor.

The “sealed” keyword means any subclass of “Dimension” must be declared in the same file as it. This way, the compiler can know, without doubt, all subclasses of Dimension, and, thus, perform checks when they get used in match statements. Leave “case Column” out of “update”, and you’ll get this warning:

MatrixInt.scala:70: warning: match is not exhaustive!
missing combination Column

what match {
^
one warning found

The usage won’t differ much either, and, in fact, won’t require any import beyond the object-companion’s contents:

scala> val a = unitMatrixInt(3)
a: teste3.MatrixInt =
/ 1 0 0 \
| 0 1 0 |
\ 0 0 1 /

scala> a(Row,2) = a.row(0)
res15: teste3.MatrixInt =
/ 1 0 0 \
| 0 1 0 |
\ 1 0 0 /

You might wonder, then, why ever use an Enumeration instead of case objects. As a matter of fact, case objects do have advantages many times, such as here. The Enumeration class, though, has many Collection methods, such as elements (iterator on Scala 2.8), which returns an Iterator, map, flatMap, filter, etc.

We have seen three different ways of implementing our “update” method for rows and columns, each with different characteristics. The very thing that made symbols a bad choice for us – the inability to define what symbols are possible in advance – can be a necessary characteristic in other situations. This goes to the very heart of Scala’s philosophy: do not dictate a solution, but let the programmer find and use the optimal one for his needs.

In the next installment we’ll look into computing cofactors, determinants, inverse matrices and then use all that to solve linear equations. I’m afraid there won’t be almost anything new, but, on the positive side, it shouldn’t take long either, and we’ll end up with a reasonable complete Matrix class.

After that, we’ll look into branching Matrix into other types, and the trade-offs we must make.

Thursday, June 25, 2009

Matrices and Fibonacci 3

Continuing our series on creating a Matrix class for Scala, now that we have the basics laid down, let’s go back to the Fibonacci problem. We need to do exponentiation of matrices. To do that, we first need to get multiplication going.

While I have no intention of explaining how to multiply matrices – the web is full of resources if you aren’t familiar with it – it’s important to restate the problem:

To find element of (row i, column j) of the matrix resulting from AxB, I need to multiply the nth element of row i of matrix A by the nth element of column j of matrix B, for all elements of said row and column, and then add the results. An initial, misguided, attempt might go like this:


(for(row <- 0 to rows
col <- 0 to B.cols
) yield A(i,col)*B(row,j)) reduceLeft (_+_)


The problem here is that you are multiplying each element of the row by every element of the column. What we want to do is to “pair” the first element of the row with the first element of the column, the second element of the row with the second of the column, and so forth, and then multiply each pair. This pairing is called “zip” in Scala. You’ll find “zip” – and the occasional variation, such as zipWithIndex – on sequence-like collections, such as List and Array. The zipping thing goes like this:


_row(i) zip other._col(j)


That will create an Array of tuples of Int: Array[Tuple2[Int,Int]], or simply Array[(Int,Int)]. Tuples can be thought of as a small “package” of objects. It’s not a collection – you do not traverse it or iterate through it. Also, each component of it has its own type, as you can see. The only thing you can do with tuples, aside from passing it around, is accessing a tuple’s element. You do this with the methods “_1”, “_2”, “_3”, etc.

Next, we need to do _1*_2 for each element of that array. You can’t do "map (_._1 * _._2)", as this translates into "map ((x1, x2) => x1._1 * x2._2)", but map only passes a single argument: the tuple. You can do this, then: "map (t => t._1 * t._2)". You can, but I won’t… Instead, I’ll pass "(_*_)", an anonymous function that takes two parameters and return the result of their multiplication, to a factory of functions.

Factories, as anyone familiar with OOP might know, are methods which create new instances of objects. So, what that has to do with a function? Well, a function in Scala is an object that implements the method “apply”. You might wonder if MatrixInt isn’t a function, then, as it implements apply itself. In fact, it could be one, if we chose to extend MatrixInt with the trait Function2 (two arguments to apply).

At any rate, we’ll use a function factory to produce a new instance of a function, based on (_*_). This factory, called Function.tupled (method “tupled” on singleton object Function), takes a function which receives multiple parameters, and converts it into a function which receives a single tuple instead. After that, we can apply the reduceLeft. Our computation, then, will go like this:


this._row(i) zip other._col(j) map Function.tupled(_*_) reduceLeft (_+_)


The method, then, is thus defined:


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


If we do away with the check to see if the matrices are compatible in size, and make allowance for the verbosity of the identifiers being used (like row and col instead of I and j), it comes down to a one-liner. Truly:


def *(o: MatrixInt) = new MatrixInt(createArrays(rows, o.cols, (i,j) => _row(i) zip o._col(j) map Function.tupled(_*_) reduceLeft (_+_)))


The only thing left to do is exponentiation. We could do it as a one-liner too:


def ^(n: Int) = (0 until n).foldLeft(unitMatrixInt(rows))((acc,_) => acc*this)


The method foldLeft is similar to reduce. But where reduce can be thought of as (((a+b)+c)+d)+e, foldLeft is more like this:


var acc = ?
acc = op(acc,a)
acc = op(acc,b)
acc = op(acc,c)
acc = op(acc,d)
acc = op(acc,e)


The most important difference between foldLeft and reduceLeft is that the result of foldLeft, and its accumulator, might be of a different type than the elements the operation is applied upon. With reduceLeft, once you applied acc = op(a,b), if acc where of a different type, then op(acc,c) would be a different method. It might work with some implicits or subtypes, but not for entirely different types.

Now, I’m not using the range for nothing, except counting the number of times I want to multiply the matrix. Acc is initialized with “unitMatrixInt(rows)”, which would be a method producing a “unit” matrix of the appropriate size. The “1” of matrices. Then I keep multiplying “this” against this accumulator.

I could create a list of copies of “this”, and then multiply them with reduceLeft too. It would go like this:


def **(n: Int) = List.fill(n)(this) reduceLeft (_*_)


Much simpler, but I wanted to mention foldLeft somehow. :-)

While simplicity has its merits, there is a simple optimization when doing exponentiation by taking advantage of common factors. I went for it in my final implementation:


def **(n: Int): MatrixInt =
if (rows != cols n < 0 =""> unitMatrixInt(rows)
case 1 => this
case 2 => this * this
case odd if odd % 2 == 1 => this ** (odd - 1) * this
case even => this ** (even / 2) ** 2
}


I'll make a note on the method name, here. I chose ** as this is one of the common names for such an operation on computer languages. Another, less common, name is ^. I prefer the later, but it has a very low precedence. In fact, this very method would be broken as written, because "(odd - 1) * this" would have precedence over "this ** (odd - 1)", which would completely invert the logic and cause a compilation error.

Now, aside from , ^, &, <>, = and !, :, +, -, * and / and % (ordered by precedence, by the way), all other special characters have a higher precedence. If I called it "@", for instance, then "a * b @ c" would behave as one would expect it too. Unfortunately, no one would expect "@" to be an exponentiation character. But if I were to make extensive use of this library in formulas and expressions, I might well go for something like that, to prevent unexpected errors due to precedence.

Anyway, we need to define unitMatrixInt too. We’ll do so in a companion object. This way, users of this package can create unix matrices if they want to. To create this method, we’ll need a more general factory method, though. We already have most of what we need in createArrays.

But before I get into that, I want to revise my class definition. Presently, this is what we have:


class MatrixInt(array: Array[Array[Int]]) {
private val self = clone(array)
}


I receive a matrix, and then replicate it internally, so I can be sure it won’t be changed outside my class – and I take care of not changing it inside. But there is a downside to it: “array” continues to be part of my object (see footnote 1)! It becomes a private reference inside my object, and this has some unfortunate consequences. For one thing, since each instance of my class makes reference to the arrays used to build them, those arrays never get deallocated. Now think about all references to “new MatrixInt” inside the class. I use that for every operator. That means I keep TWO copies of a MatrixInt’s array in memory, even when the copy used for construction gets immediately discarded.

This presents a conundrum.On one hand, if the instances of my class are immutable, then I need the full array to begin with. On the other hand, if I receive an array from an external source, I need to copy it before using, as the original array might get changed by whoever passed it to me.

We solve this with the help of a private constructor and an object companion. A private constructor is, of course, a constructor that cannot be accessed from outside the class. The private constructor thing is easy:


class MatrixInt private (array: Array[Array[Int]]) {


Now, how to access it? We cannot define another constructor inside MatrixInt with the same signature – receiving an array of array of Int as parameter – because it would conflict with the first. What we do, then, is use an object companion.

Simply put, an object companion is a singleton object – declared with the “object” keyword – which has the same name as a class. That grants both the object and the class visibility of each other’s private definitions.

We’ll then define an “apply” method in that object companion. As you may remember, the “apply” method receives special treatment from the Scala compiler, which enables you to do away with the name – apply – and just pass the parameters to the instance of your object. Since we have to call our object companion “MatrixInt”, then MatrixInt(…) will be converted by Scala into MatrixInt.apply(…). It looks like the instantiation of a new object, only it lacks the “new” keyword.

Let’s define it, then. First, let’s remove the “val self” from the class, and just receive it in its private constructor. As that constructor will only be accessed through MatrixInt’s object companion, that will be safe.


class MatrixInt private (self: Array[Array[Int]]) {
import MatrixInt._


The import statement makes it possible to access the methods of the object companion without preceding them with “MatrixInt.”. Being a companion object makes the private methods accessible, but does not bring them into scope automatically.

Now, for the object companion, let’s create it, the factory for MatrixInt, and a couple of helpful factories. We’ll also transfer the “clone” and “createArrays” methods to it, as this method does not make reference to a MatrixInt’s own data, so it doesn’t need to be part of each instance. The definition of toArray needs “MatrixInt.” prepended to “clone” to avoid an ambiguous reference, but nothing else should change.


object MatrixInt {
import MatrixInt._
implicit def toMatrixInt(m : Array[Array[Int]]) = MatrixInt(m)

def apply(array: Array[Array[Int]]): MatrixInt = new MatrixInt(clone(array))
def apply(rows: Int, cols: Int): MatrixInt = MatrixInt(rows, cols, 0)
def apply(rows: Int, cols: Int, value: Int): MatrixInt =
new MatrixInt(createArrays(rows, cols, ((_,_) => value)))
def apply(rows: Int, cols: Int, f: (Int,Int) => Int): MatrixInt =
new MatrixInt(createArrays(rows, cols, f))

def unitMatrixInt(n: Int) = MatrixInt(n, n, (row,col) => (if (row == col) 1 else 0))

protected def createArrays(rows: Int, cols: Int, f: (Int, Int) => Int) =
for((i: Int) <- (0 until rows).toArray) yield for((j: Int) <- (0 until cols).toArray) yield f(i,j) protected def clone(a: Array[Array[Int]]) = createArrays(a.size, a(0).size, (row, col) => a(row)(col))
}


We can now create a MatrixInt by using MatrixInt(Array(Array(…))), but we didn’t – and won’t – change that in the class. The companion object factory produces a clone of the array, just as we wanted. But the arrays produced by MatrixInt’s methods don’t get exposed, so there is no risk in using them. So, by keep using the private constructor inside MatrixInt itself, we prevent unnecessary copies of the arrays produced by each of MatrixInt’s operators.

We’ll finish today with Fibonacci, but we’ll continue later on. I’ll add some methods for completeness to MatrixInt, and then look into MatrixDouble and solve some linear equations with it.

After that, we’ll briefly explore ways to make our Matrix class parametrizable, though I’m still in doubt about performance. Still, I’m reading the comments, and I plan to address them one way or another in these posts.

So, the Fibonacci!


def fibonacci(n: Int) = (MatrixInt(Array(Array(1,1),Array(1,0))) ** n)(1,0)


Footnote 1:

A reader has brought to my attention that, actually, "array" is not being preserved. I haven't found in Scala specification anything about this, so I'm putting it down to optimization. But because it is not explicitly documented, I advise against taking advantage of it.

Do note, please, that if you declare "array" to be protected, val or var, or if any method makes use of it, then it will be preserved in the object.

Wednesday, June 24, 2009

Matrices 2

I’ll continue today with a little series for beginners. On my first post I forgot something important. One of the very first things I do with a class is create a toString method. This way, I get instant feedback on how my class is doing, as I paste it into the REPL, and try out the methods I defined.

So, let’s think of a toString method. Since all data in my class is stored in the private val “self”, we could simply delegate the task of producing a string:


override def toString = self.toString


That might do at first, but it would be unreasonable to expect something better formatted than that. First, let’s see if we can’t get the output in multiple lines, using mkString:


override def toString = self mkString “\n”


Unfortunately, the output of that looks like this:


res6: String =
[I@196a753
[I@1c36f46


This wouldn’t happen with a List, but an Array is implemented to be compatible with Java, and Java had Arrays have some intrinsic limitations, as they were present in the language from the very beginning, much before generics, and some design decisions about them proved to be unfortunate. Anyway, let’s fix that:


override def toString = self map (_ mkString “ “) mkString “\n”


Better, but we can improve on that still. For one thing, it would be much better if the columns were aligned. Now, to align the columns, we first must know how big a column can get. Or, in other words, we must know the size of the number with the biggest string representation in the matrix.

So, summing it up, we have to look through all of the matrix elements, get their sizes, and choose the biggest one. Getting a number’s size is easy: number.toString.size. Now, if we had a list, this would be easy, but we have an array of arrays. When we have nested collections and we want to unnest them, we use flatMap. For instance:


scala> Array(Array(15,0),Array(2,-30)) flatMap (x => x)
res19: Array[Int] = Array(15, 0, 2, -30)


The flatMap operation we used, (x => x), simply returns what it was given. If all you want is unnest the collections, this suffices. But, actually, we want the number sizes. Since “x” is an Array, not a number, we can’t do flatMap (x => x.toString.size). But we can replace the numbers with their sizes, and then unnest the arrays, like this:


scala> Array(Array(15,0),Array(2,-30)) flatMap (_ map (_.toString.size))
res20: Array[Int] = Array(2, 1, 1, 3)


Next, we want to get the biggest one. We have a collection of number sizes, and we want to reduce it to the single biggest one. We have two such methods in Scala: reduceLeft and reduceRight. Since “max” is commutative, both will produce the same result. But reduceLeft has a better performance, as it iterates over the collection in the order given. The final code should look like this, then:


val maxNumSize = self flatMap (_ map (_.toString.size)) reduceLeft (_ max _)


Pretty easy, eh? My first version had three lines, and was much more confusing to read. It helps stating what your problem is – you’ll often see you can do exactly that in Scala. So, we have a size, what do we do with it? Well, let’s just make sure all numbers are printed to that size. If we format them with “%2d”, for instance, every number will take at least two characters, with the left side being filled with a space as needed. To produce such formatted numbers, we’ll make use of RichString’s “format” method:


override def toString = {
val maxNumSize = self.projection flatMap (_ map (_.toString.size)) reduceLeft (_ max _)
val numFormat = "%"+maxNumSize+"d"
self map (_ map (numFormat format _) mkString “ “) mkString “\n”
}


That’s it! I added a “projection” there to make it a bit more performatic, as it doesn’t produce intermediate arrays. On Scala 2.8, that method gets replaced with “view”, with “projection” remaining as deprecated. Well, that was so easy I want to do more. First, I want to produce little vertical bars on the sides, to make clear the matrix boundaries, in case we are printing them side by side. That just requires changing the inner mkString to mkString(“| “, “ “, “ |”). But I also want to bend the bars up and down, to indicate the top and bottom lines, as well as produce something special for one-line matrices.

For all of them, we are just making small changes in the inner mkString. We have one for the top line, one for the intermediate ones, and one for the bottom line. We could make a list of all these lines, and then use mkString on the list, but we have to be careful, because the “bottom” line and the “intermediate” lines might not exist, and the top line has two different formats depending on whether the matrix has more than one line or not. Taking that into account, here is a much larger, but, hopefully, nicer toString method:


override def toString = {
val maxNumSize = self.projection flatMap (_ map (_.toString.size)) reduceLeft (_ max _)
val numFormat = "%"+maxNumSize+"d"
def formatNumber(n: Int) = numFormat format n
val top = rows match {
case 1 => _row(0) map formatNumber mkString ("< ", " ", " >")
case _ => _row(0) map formatNumber mkString ("/ ", " ", " \\")
}
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"
}


In a parting note, someone suggested one technique to be able to parametrize the class. There are many, but they are slower than using "AnyVal" types directly. If I happen to be mistaken, please do post about it at length!

Tuesday, June 23, 2009

Matrices and Fibonacci 1

When I was beginning to learn Scala, I decided one of my projects would be rewrite this stuff in Scala. The idea behind that is to compute Fibonacci Numbers using matrices. I won’t bother describing how it is done, as this is well covered on both of those links. Suffice to say that I need to computer the nth power of a 2x2 matrix of integers. This project was all but forgotten for a long time…

So it happened that just yesterday I was cleaning my to-do list, and I came upon that. I took a second look into it, and decided to take the plunge. Now, Raganwald uses an optimization which enables him to do the operations on vectors, but, by now, I actually thought the whole exercise was trivial and pointless.

I took a brief look into matrix libraries for Scala. I found the Scalala (as in Scala Linear Algebra) library, which does everything I need, very efficiently at that. My problem would be reduced to something like this (untested):


def fib(n) = n match {
case 0 => throw new IllegalArgumentException
case 1 => 1
case _ =>
val m = DenseMatrix(2,2)(1); m(1,1) = 0
val mPowered = m :^ n
mPowered(0,1)
}


Boring, wasn’t it? Yeah, I thought so too. Instead, I decided to write my own Matrix library. Pretty trivial stuff, but at least I’d take something out of it, I thought. So, I’ll talk about that: writing my Matrix library. If that will bore you, you might as well stop right here. I’ll show, though, some of the design decisions and constrains I faced. For someone new to Scala, it might be interesting read.

Well, then, where to start? With a class, of course. Ideally, I’d like to write something like this:


class Matrix[T]


This way, I could use matrices of Double for most use cases, of Ints or Longs for things like the Fibonacci problem, and BigDecimal for tough stuff. Unfortunately, we can’t do that. We need access to operators such as +, - and * over type T, and that means I need an upper bound for T defining these operators. Unfortunately, these types do not share an interface defining those elements. You see, they are not really “classes” as far as Java is concerned, and not being a proper class, they can’t share an abstract class or interface with these methods. Well, BigDecimal is a class, alright, but that’s it.

Because of this, I decided, instead, to define my class as MatrixInt. This was enough for Fibonacci, but later I created a MatrixDouble as well, just to finish implementing the most basic methods such a class would need.

There is a way around that. I can get away with a structural type, but that would be so slow as to defeat the whole purpose of a Matrix class. So I’ll skip it altogether.

The next design decision is how to store my data, and what constructor to use. Generally speaking, there are two strategies for storing a matrix: as a dense matrix or as a sparse matrix. A dense matrix is simply stored as an addressable, fixed-size, contiguously allocated buffer. Or, in other words, an array. A sparse matrix is used for matrices that are very large in their dimensions, but are mostly populated with zeros. Performance can go either way, depending on what kind of data and algorithms you are using those matrices for. For Fibonacci, dense is the way to go. Since it is also the easier way to go, I chose that way.

We also need to know how many dimensions the matrices may have. For simplicity sake, I chose bi-dimensional matrices, and that’s it. With all that in mind, my class thus began:


class MatrixInt(array: Array[Array[Int]]) {


I decided to take an array of proper dimensions as my constructor, instead of providing something special myself. With Scala 2.8, you can pass Array.ofDim(2,2) and be done with it. We’ll see something similar further down too.

Now, one important thing to consider is that an Array is a mutable object. I have decided to make my matrices immutable, so that I could use them just like I use Ints or Doubles. But, to ensure that is true, I can’t really on the array passed with the constructor, as it might change. So, first thing I do is copy it. Here:


private def clone(a: Array[Array[Int]]) =
for((i: Int) <- (0 until a.size).toArray)
yield for((j: Int) <- (0 until a(0).size).toArray)
yield a(i)(j)

private val self = clone(array)


This clone array method will help in places other than the constructor. For instance, if the user wants the Array representation of the matrix, I have to clone my internal representation before passing it on. But there are more interesting things about it. Let’s consider it…

One thing to notice about it, is that instead of having a single for/yield statement, I have two, the second being fed to the first’s yield. It is built this way because I want to create an Array of Arrays. Each yield will result in one collection. Two collections, two yield statements.

Then we have ranges: (0 until a.size) and (0 until a(0).size). A statement like (0 until 5) will produce a range which, when iterated over, will yield 0, 1, 2, 3 and 4. It can do a few other things too, but let’s not be bothered by them. Note that if you want an “inclusive”, range, one which will include the upper limit (5, in the example), you use “to” instead of “until”. Or you create a new range from the existing one with “inclusive”.

Next, “.toArray”. Though a range will do as far as generating indices, I want yield to return an Array. If you want yield to return collection X, pass a collection X to the for statement. Therefore, “.toArray” to convert the range into an array.

The final thing to consider is the last yield statement. Yield a(i)(j). As I started writing my methods, I got very similar patterns, with only the last yield changing. So I started to consider how to reuse that code. The thing to consider, here, is that the intent of this double for/yield pattern is to create an array of arrays of the proper size, and populate it with a function of my choosing. Thinking about it that way, a refactoring might be done as follow:


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

protected def clone(a: Array[Array[Int]]) =
createArrays(a.size, a(0).size, (row, col) => a(row)(col))


Here I changed from private to protected, as I decided if I ever subclassed it, there would be no point in hiding these methods. I do not publish them, though, because they are helper methods, instead of methods acting upon the object itself.

Next, a.size and a(0).size were a bit cryptic, weren’t they? When we work on the matrices algorithms, it would be preferable to call them by what they are:


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


We would also like to be able to get a particular row or a particular column. It seems getting a particular row would be easy. Just return the subarray:


def row(n : Int): Array[Int] = self(n)


There is a serious problem with that, though. I’m returning a mutable object from my own internal data! We could solve that with clone, but clone only worked for bi-dimensional arrays. Well, before we create new clone, let’s consider a simpler alternative:


private def _row(n: Int): Array[Int] = self(n)
def row(n: Int): Array[Int] = _row(n) map (x => x)


The private definition is here so we can avoid creating new arrays unnecessarily. Internally, we’ll make sure we never change the value of our arrays. The public method creates a new array by the magic of “map”. The method “map” is, in fact, used when you do a for/yield, along with flatMap when the for/yield has multiple iterators. So, in fact, the createArrays method is actually implemented as using map.

A map will take a collection and produce a new collection with a given mapping. In this case, our mapping is (x => x), or, in other words, given x, return x, producing a copy of that array. We define cols in a likewise fashion:


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


In this case, the private definition and the public definition are the same. Still, I decided to keep a private definition so that I could use _row and _col throughout my class. We’ll also define a method to access a particular element of our matrix:


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


The method “apply” has a special meaning in Scala, as the compiler translates object(parms) into object.apply(parms). This means that I will be able to access a particular element of a matrix m with “m(i,j)”. Of course, because a matrix has bounds, I must test for them, and throw an appropriate exception if need be.

Note that I do not declare an @throws annotation, though I might if I wished to. Also, note that I’m returning the result of the if/else expression, and its type should be Int. But in the if statement, I’m throwing an exception, so what gives?

As it happens, throwing an exception has type Nothing, which is subtype to everything. I’m not sure if the compiler actually assigns a type to “throw”, as this is a keyword. But if I used Scala’s library “error” function instead, the question of the type would be important. For methods and functions which never return, such as error, declaring its type to be Nothing makes it possible to make use of them inside statements whose return type is important.

We’ll stop today with just three more definitions. Let’s multiply a matrix by an integer, producing a new matrix, as well as add and subtract two arrays. Here is how:


def *(n: Int): MatrixInt =
new MatrixInt(createArrays(rows, cols, (row, col) => this(row,col) * n))

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

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