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")
}
Posts
No comments:
Post a Comment