Merge pull request #154 from kylecorry31/matrix-rework

Matrix rework

Author: kylecorry31

src/main/kotlin/com/kylecorry/sol/math/algebra/LinearAlgebra.kt

@@ -1,7 +1,6 @@
 package com.kylecorry.sol.math.algebra
 
 import com.kylecorry.sol.math.SolMath
-import com.kylecorry.sol.shared.ArrayUtils.swap
 import kotlin.math.absoluteValue
 import kotlin.math.min
 import kotlin.math.sqrt
@@ -13,7 +12,7 @@ object LinearAlgebra {
             throw Exception("Matrix 1 columns must be the same size as matrix 2 rows")
         }
 
-        val product = createMatrix(mat1.rows(), mat2.columns()) { _, _ -> 0f }
+        val product = Matrix.zeros(mat1.rows(), mat2.columns())
         for (r in 0 until mat1.rows()) {
             for (otherC in 0 until mat2.columns()) {
                 var sum = 0.0f
@@ -36,7 +35,7 @@ object LinearAlgebra {
             throw Exception("Matrix 1 rows must be the same size as matrix 2 rows")
         }
 
-        return createMatrix(mat1.rows(), mat1.columns()) { row, col ->
+        return Matrix.create(mat1.rows(), mat1.columns()) { row, col ->
             mat1[row, col] - mat2[min(row, mat2.rows() - 1), min(col, mat2.columns() - 1)]
         }
     }
@@ -54,13 +53,13 @@ object LinearAlgebra {
             throw Exception("Matrix 1 rows must be the same size as matrix 2 rows")
         }
 
-        return createMatrix(mat1.rows(), mat1.columns()) { row, col ->
+        return Matrix.create(mat1.rows(), mat1.columns()) { row, col ->
             mat1[row, col] + mat2[min(row, mat2.rows() - 1), min(col, mat2.columns() - 1)]
         }
     }
 
     fun add(mat1: Matrix, value: Float): Matrix {
-        return createMatrix(mat1.rows(), mat1.columns()) { row, col ->
+        return Matrix.create(mat1.rows(), mat1.columns()) { row, col ->
             mat1[row, col] + value
         }
     }
@@ -74,13 +73,13 @@ object LinearAlgebra {
             throw Exception("Matrix 1 rows must be the same size as matrix 2 rows")
         }
 
-        return createMatrix(mat1.rows(), mat1.columns()) { row, col ->
+        return Matrix.create(mat1.rows(), mat1.columns()) { row, col ->
             mat1[row, col] * mat2[min(row, mat2.rows() - 1), min(col, mat2.columns() - 1)]
         }
     }
 
     fun multiply(mat1: Matrix, scale: Float): Matrix {
-        return createMatrix(mat1.rows(), mat1.columns()) { row, col ->
+        return Matrix.create(mat1.rows(), mat1.columns()) { row, col ->
             mat1[row, col] * scale
         }
     }
@@ -94,7 +93,7 @@ object LinearAlgebra {
             throw Exception("Matrix 1 rows must be the same size as matrix 2 rows")
         }
 
-        return createMatrix(mat1.rows(), mat1.columns()) { row, col ->
+        return Matrix.create(mat1.rows(), mat1.columns()) { row, col ->
             mat1[row, col] / mat2[min(row, mat2.rows() - 1), min(col, mat2.columns() - 1)]
         }
     }
@@ -104,50 +103,60 @@ object LinearAlgebra {
     }
 
     fun transpose(mat: Matrix): Matrix {
-        return createMatrix(mat.columns(), mat.rows()) { row, col ->
+        return Matrix.create(mat.columns(), mat.rows()) { row, col ->
             mat[col, row]
         }
     }
 
     fun map(mat: Matrix, fn: (value: Float) -> Float): Matrix {
-        return createMatrix(mat.rows(), mat.columns()) { row, col ->
+        return Matrix.create(mat.rows(), mat.columns()) { row, col ->
             fn(mat[row, col])
         }
     }
 
     fun mapRows(mat: Matrix, fn: (row: FloatArray) -> FloatArray): Matrix {
-        return mat.map { fn(it.toFloatArray()).toTypedArray() }.toTypedArray()
+        val copy = mat.clone()
+        val temp = FloatArray(mat.columns())
+        for (row in 0 until mat.rows()) {
+            mat.getRow(row, temp)
+            copy.setRow(row, fn(temp))
+        }
+        return copy
     }
 
     fun mapColumns(mat: Matrix, fn: (row: FloatArray) -> FloatArray): Matrix {
         return mapRows(mat.transpose(), fn).transpose()
     }
 
     fun sum(mat: Matrix): Float {
-        return mat.sumOf { it.sum().toDouble() }.toFloat()
+        return mat.sum()
     }
 
     fun sumColumns(mat: Matrix): Matrix {
         return sumRows(mat.transpose()).transpose()
     }
 
     fun sumRows(mat: Matrix): Matrix {
-        return createMatrix(mat.rows(), 1) { row, _ ->
-            mat[row].sum()
+        val temp = FloatArray(mat.columns())
+        return Matrix.create(mat.rows(), 1) { row, _ ->
+            mat.getRow(row, temp)
+            temp.sum()
         }
     }
 
     fun max(mat: Matrix): Float {
-        return mat.maxOf { it.max() }
+        return mat.max()
     }
 
     fun maxColumns(mat: Matrix): Matrix {
         return maxRows(mat.transpose()).transpose()
     }
 
     fun maxRows(mat: Matrix): Matrix {
-        return createMatrix(mat.rows(), 1) { row, _ ->
-            mat[row].max()
+        val temp = FloatArray(mat.columns())
+        return Matrix.create(mat.rows(), 1) { row, _ ->
+            mat.getRow(row, temp)
+            temp.max()
         }
     }
 
@@ -159,7 +168,7 @@ object LinearAlgebra {
         val det = determinant(m)
         if (SolMath.isZero(det)) {
             // No inverse exists
-            return createMatrix(m.rows(), m.columns(), 0f)
+            return Matrix.zeros(m.rows(), m.columns())
         }
         return adjugate(m).transpose().divide(determinant(m))
     }
@@ -171,7 +180,7 @@ object LinearAlgebra {
 
         var colMultiplier: Int
         var rowMultiplier: Int
-        return createMatrix(m.rows(), m.columns()) { r, c ->
+        return Matrix.create(m.rows(), m.columns()) { r, c ->
             rowMultiplier = if (r % 2 == 0) {
                 1
             } else {
@@ -208,7 +217,7 @@ object LinearAlgebra {
     }
 
     fun cofactor(m: Matrix, r: Int, c: Int): Matrix {
-        return createMatrix(m.rows() - 1, m.columns() - 1) { r1, c1 ->
+        return Matrix.create(m.rows() - 1, m.columns() - 1) { r1, c1 ->
             val sr = if (r1 < r) {
                 r1
             } else {
@@ -227,20 +236,20 @@ object LinearAlgebra {
         val rows = m.rows()
         val cols = m.columns()
 
-        val q = createMatrix(rows, cols) { _, _ -> 0f }
-        val r = createMatrix(cols, cols) { _, _ -> 0f }
+        val q = Matrix.zeros(rows, cols)
+        val r = Matrix.zeros(cols, cols)
 
         for (j in 0 until cols) {
-            var v = m.column(j)
+            var v = Matrix.row(values = m.getColumn(j))
 
             for (i in 0 until j) {
-                val qi = q.column(i)
+                val qi = Matrix.row(values = q.getColumn(i))
                 r[i, j] = qi.dot(v.transpose())[0, 0]
                 v = v.subtract(qi.multiply(r[i, j]))
             }
 
             r[j, j] = norm(v)
-            q.setColumn(j, v.divide(r[j, j])[0])
+            q.setColumn(j, v.divide(r[j, j]).getRow(0))
         }
 
         return q to r
@@ -250,12 +259,12 @@ object LinearAlgebra {
      * Returns a column matrix of the diagonal of the matrix
      */
     fun diagonal(m: Matrix): Matrix {
-        return createMatrix(1, min(m.rows(), m.columns())) { i, _ ->
+        return Matrix.create(1, min(m.rows(), m.columns())) { i, _ ->
             m[i, i]
         }
     }
 
-    fun eigenvalues(m: Matrix, tolerance: Float = 1e-12f, maxIterations: Int = 1000): Array<Float> {
+    fun eigenvalues(m: Matrix, tolerance: Float = 1e-12f, maxIterations: Int = 1000): FloatArray {
         var old = m
         var new = m
 
@@ -269,25 +278,19 @@ object LinearAlgebra {
             iterations++
         }
 
-        return diagonal(new)[0]
+        return diagonal(new).getRow(0)
     }
 
     fun norm(m: Matrix): Float {
-        return sqrt(
-            m.sumOf { values ->
-                values.sumOf { value ->
-                    value * value.toDouble()
-                }
-            }
-        ).toFloat()
+        return m.norm()
     }
 
-    fun solveLinear(a: Matrix, b: Array<Float>): Array<Float> {
+    fun solveLinear(a: Matrix, b: FloatArray): FloatArray {
         require(a.rows() == a.columns()) { "Matrix must be square" }
         require(a.rows() == b.size) { "Matrix rows must be the same size as the vector" }
 
         val n = a.columns()
-        val ab = a.appendColumn(b.toFloatArray())
+        val ab = a.appendColumn(b)
 
         // Convert to row echelon form
         for (i in 0 until n) {
@@ -298,7 +301,7 @@ object LinearAlgebra {
                 }
             }
 
-            ab.swap(i, maxRow)
+            ab.swapRows(i, maxRow)
 
             for (j in i + 1 until n) {
                 val factor = ab[j, i] / ab[i, i]
@@ -317,13 +320,13 @@ object LinearAlgebra {
             }
         }
 
-        return x.toTypedArray()
+        return x
     }
 
-    fun leastNorm(a: Matrix, b: Array<Float>): Array<Float> {
+    fun leastNorm(a: Matrix, b: Array<Float>): FloatArray {
         val (q, r) = qr(a.transpose())
-        val y = q.dot(r.inverse().transpose()).dot(createMatrix(b.size, 1) { i, _ -> b[i] })
-        return y.transpose()[0]
+        val y = q.dot(r.inverse().transpose()).dot(Matrix.create(b.size, 1) { i, _ -> b[i] })
+        return y.getColumn(0)
     }
 
     /**
@@ -336,17 +339,17 @@ object LinearAlgebra {
     fun leastSquares(a: Matrix, b: Vector): Vector {
         val isUnderdetermined = a.rows() < a.columns()
         if (isUnderdetermined) {
-            return leastNorm(a, b)
+            return leastNorm(a, b).toTypedArray()
         }
 
         val jt = a.transpose()
         val jtj = jt.dot(a)
         val jtr = jt.dot(b.toColumnMatrix())
-        return solveLinear(jtj, jtr.transpose()[0])
+        return solveLinear(jtj, jtr.getColumn(0)).toTypedArray()
     }
 
     fun appendColumn(m: Matrix, col: FloatArray): Matrix {
-        return createMatrix(m.rows(), m.columns() + 1) { r, c ->
+        return Matrix.create(m.rows(), m.columns() + 1) { r, c ->
             if (c < m.columns()) {
                 m[r, c]
             } else {
@@ -356,7 +359,7 @@ object LinearAlgebra {
     }
 
     fun appendColumn(m: Matrix, value: Float): Matrix {
-        return createMatrix(m.rows(), m.columns() + 1) { r, c ->
+        return Matrix.create(m.rows(), m.columns() + 1) { r, c ->
             if (c < m.columns()) {
                 m[r, c]
             } else {
@@ -366,7 +369,7 @@ object LinearAlgebra {
     }
 
     fun appendRow(m: Matrix, row: FloatArray): Matrix {
-        return createMatrix(m.rows() + 1, m.columns()) { r, c ->
+        return Matrix.create(m.rows() + 1, m.columns()) { r, c ->
             if (r < m.rows()) {
                 m[r, c]
             } else {
@@ -376,7 +379,7 @@ object LinearAlgebra {
     }
 
     fun appendRow(m: Matrix, value: Float): Matrix {
-        return createMatrix(m.rows() + 1, m.columns()) { r, c ->
+        return Matrix.create(m.rows() + 1, m.columns()) { r, c ->
             if (r < m.rows()) {
                 m[r, c]
             } else {