## 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))`

`  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> mres5: Array[Array[Int]] = Array(Array(5, -4), Array(6, -5))scala> rres6: Array[Array[Int]] = Array(Array(2), Array(1))scala> new LinearEquations(m, r)res2: LinearEquations =5x(0) + -4x(1) = 26x(0) + -5x(1) = 1scala> res2.solveInverseres3: java.lang.String =5*6 + -4*7 = 26*6 + -5*7 = 1scala> res2.solveCramerres4: java.lang.String =5*6 + -4*7 = 26*6 + -5*7 = 1`

And now, the source code:
`package com.blogspot.dcsobral.matrixclass 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")}`