WIP Matrix refactor

examples/src/main/kotlin/space/kscience/kmath/structures/typeSafeDimensions.kt

@@ -5,7 +5,7 @@ import space.kscience.kmath.dimensions.D3
 import space.kscience.kmath.dimensions.DMatrixContext
 import space.kscience.kmath.dimensions.Dimension
 
-private fun DMatrixContext<Double>.simple() {
+private fun DMatrixContext<Double, *>.simple() {
     val m1 = produce<D2, D3> { i, j -> (i + j).toDouble() }
     val m2 = produce<D3, D2> { i, j -> (i + j).toDouble() }
 
@@ -17,7 +17,7 @@ private object D5 : Dimension {
     override val dim: UInt = 5u
 }
 
-private fun DMatrixContext<Double>.custom() {
+private fun DMatrixContext<Double, *>.custom() {
     val m1 = produce<D2, D5> { i, j -> (i + j).toDouble() }
     val m2 = produce<D5, D2> { i, j -> (i - j).toDouble() }
     val m3 = produce<D2, D2> { i, j -> (i - j).toDouble() }

kmath-commons/src/main/kotlin/space/kscience/kmath/commons/linear/CMMatrix.kt

@@ -99,8 +99,8 @@ public object CMLinearSpace : LinearSpace<Double, RealField> {
         ArrayRealVector(array).wrap()
     }
 
-    private fun RealMatrix.wrap(): CMMatrix = CMMatrix(this)
+    internal fun RealMatrix.wrap(): CMMatrix = CMMatrix(this)
-    private fun RealVector.wrap(): CMVector = CMVector(this)
+    internal fun RealVector.wrap(): CMVector = CMVector(this)
 
     override fun buildVector(size: Int, initializer: RealField.(Int) -> Double): Vector<Double> =
         ArrayRealVector(DoubleArray(size) { RealField.initializer(it) }).wrap()

kmath-commons/src/main/kotlin/space/kscience/kmath/commons/linear/CMSolver.kt

@@ -27,7 +27,7 @@ public fun CMLinearSpace.solve(
     a: Matrix<Double>,
     b: Matrix<Double>,
     decomposition: CMDecomposition = CMDecomposition.LUP
-): CMMatrix = solver(a, decomposition).solve(b.toCM().origin).asMatrix()
+): CMMatrix = solver(a, decomposition).solve(b.toCM().origin).wrap()
 
 public fun CMLinearSpace.solve(
     a: Matrix<Double>,
@@ -38,4 +38,4 @@ public fun CMLinearSpace.solve(
 public fun CMLinearSpace.inverse(
     a: Matrix<Double>,
     decomposition: CMDecomposition = CMDecomposition.LUP
-): CMMatrix = solver(a, decomposition).inverse.asMatrix()
+): CMMatrix = solver(a, decomposition).inverse.wrap()

kmath-core/src/commonMain/kotlin/space/kscience/kmath/linear/MatrixBuilder.kt

@@ -3,11 +3,30 @@ package space.kscience.kmath.linear
 import space.kscience.kmath.misc.UnstableKMathAPI
 import space.kscience.kmath.operations.Ring
 
+public class MatrixBuilder<T : Any, A : Ring<T>>(
+    public val linearSpace: LinearSpace<T, A>,
+    public val rows: Int,
+    public val columns: Int,
+) {
+    public operator fun invoke(vararg elements: T): Matrix<T> {
+        require(rows * columns == elements.size) { "The number of elements ${elements.size} is not equal $rows * $columns" }
+        return linearSpace.buildMatrix(rows, columns) { i, j -> elements[i * columns + j] }
+    }
+
+    //TODO add specific matrix builder functions like diagonal, etc
+}
 
 @UnstableKMathAPI
-public fun <T : Any> LinearSpace<T, Ring<T>>.matrix(rows: Int, columns: Int, vararg elements: T): Matrix<T> {
+public fun <T : Any, A : Ring<T>> LinearSpace<T, A>.matrix(rows: Int, columns: Int): MatrixBuilder<T, A> =
-    require(rows * columns == elements.size) { "The number of elements ${elements.size} is not equal $rows * $columns" }
+    MatrixBuilder(this, rows, columns)
-    return buildMatrix(rows, columns) { i, j -> elements[i * columns + j] }
+
+/**
+ * Build a square matrix from given elements.
+ */
+@UnstableKMathAPI
+public fun <T : Any> LinearSpace<T, Ring<T>>.square(vararg elements: T): Matrix<T> {
+    val size: Int = kotlin.math.sqrt(elements.size.toDouble()).toInt()
+    return matrix(size,size)(*elements)
 }
 
 @UnstableKMathAPI

kmath-core/src/commonTest/kotlin/space/kscience/kmath/linear/MatrixTest.kt

@@ -38,7 +38,7 @@ class MatrixTest {
         infix fun Matrix<Double>.pow(power: Int): Matrix<Double> {
             var res = this
             repeat(power - 1) {
-                res = RealLinearSpace.invoke { res dot this@pow }
+                res = LinearSpace.real.run { res dot this@pow }
             }
             return res
         }

kmath-core/src/commonTest/kotlin/space/kscience/kmath/linear/RealLUSolverTest.kt

@@ -14,12 +14,12 @@ class RealLUSolverTest {
 
     @Test
     fun testDecomposition() {
-        val matrix = Matrix.square(
-            3.0, 1.0,
-            1.0, 3.0
-        )
-
         LinearSpace.real.run {
+            val matrix = square(
+                3.0, 1.0,
+                1.0, 3.0
+            )
+
             val lup = lup(matrix)
 
             //Check determinant
@@ -31,14 +31,14 @@ class RealLUSolverTest {
 
     @Test
     fun testInvert() {
-        val matrix = Matrix.square(
+        val matrix = LinearSpace.real.square(
             3.0, 1.0,
             1.0, 3.0
         )
 
         val inverted = LinearSpace.real.inverseWithLup(matrix)
 
-        val expected = Matrix.square(
+        val expected = LinearSpace.real.square(
             0.375, -0.125,
             -0.125, 0.375
         )

kmath-core/src/commonTest/kotlin/space/kscience/kmath/structures/NumberNDFieldTest.kt

@@ -1,5 +1,6 @@
 package space.kscience.kmath.structures
 
+import space.kscience.kmath.linear.LinearSpace
 import space.kscience.kmath.nd.*
 import space.kscience.kmath.operations.Norm
 import space.kscience.kmath.operations.invoke
@@ -33,7 +34,9 @@ class NumberNDFieldTest {
     @Test
     fun testGeneration() {
 
-        val array = Structure2D.real(3, 3) { i, j -> (i * 10 + j).toDouble() }
+        val array = LinearSpace.real.buildMatrix(3, 3) { i, j ->
+            (i * 10 + j).toDouble()
+        }
 
         for (i in 0..2) {
             for (j in 0..2) {

kmath-for-real/src/commonMain/kotlin/space/kscience/kmath/real/dot.kt

@@ -9,23 +9,4 @@ import space.kscience.kmath.nd.Matrix
  */
 public infix fun Matrix<Double>.dot(other: Matrix<Double>): Matrix<Double> = LinearSpace.real.run{
     this@dot dot other
-//    require(colNum == other.rowNum) { "Matrix dot operation dimension mismatch: ($rowNum, $colNum) x (${other.rowNum}, ${other.colNum})" }
-//    val resultArray = DoubleArray(this.rowNum * other.colNum)
-//
-//    //convert to array to insure there is no memory indirection
-//    fun Buffer<Double>.unsafeArray() = if (this is RealBuffer)
-//        this.array
-//    else
-//        DoubleArray(size) { get(it) }
-//
-//    val a = this.buffer.unsafeArray()
-//    val b = other.buffer.unsafeArray()
-//
-//    for (i in (0 until rowNum))
-//        for (j in (0 until other.colNum))
-//            for (k in (0 until colNum))
-//                resultArray[i * other.colNum + j] += a[i * colNum + k] * b[k * other.colNum + j]
-//
-//    val buffer = RealBuffer(resultArray)
-//    return BufferMatrix(rowNum, other.colNum, buffer)
 }

kmath-for-real/src/commonTest/kotlin/kaceince/kmath/real/RealMatrixTest.kt

@@ -1,7 +1,7 @@
 package kaceince.kmath.real
 
-import space.kscience.kmath.linear.Matrix
+import space.kscience.kmath.linear.LinearSpace
-import space.kscience.kmath.linear.build
+import space.kscience.kmath.linear.matrix
 import space.kscience.kmath.misc.UnstableKMathAPI
 import space.kscience.kmath.real.*
 import space.kscience.kmath.structures.contentEquals
@@ -30,11 +30,11 @@ internal class RealMatrixTest {
 
     @Test
     fun testRepeatStackVertical() {
-        val matrix1 = Matrix.build(2, 3)(
+        val matrix1 = LinearSpace.real.matrix(2, 3)(
             1.0, 0.0, 0.0,
             0.0, 1.0, 2.0
         )
-        val matrix2 = Matrix.build(6, 3)(
+        val matrix2 = LinearSpace.real.matrix(6, 3)(
             1.0, 0.0, 0.0,
             0.0, 1.0, 2.0,
             1.0, 0.0, 0.0,
@@ -47,12 +47,12 @@ internal class RealMatrixTest {
 
     @Test
     fun testMatrixAndDouble() {
-        val matrix1 = Matrix.build(2, 3)(
+        val matrix1 = LinearSpace.real.matrix(2, 3)(
             1.0, 0.0, 3.0,
             4.0, 6.0, 2.0
         )
         val matrix2 = (matrix1 * 2.5 + 1.0 - 2.0) / 2.0
-        val expectedResult = Matrix.build(2, 3)(
+        val expectedResult = LinearSpace.real.matrix(2, 3)(
             0.75, -0.5, 3.25,
             4.5, 7.0, 2.0
         )
@@ -61,13 +61,13 @@ internal class RealMatrixTest {
 
     @Test
     fun testDoubleAndMatrix() {
-        val matrix1 = Matrix.build(2, 3)(
+        val matrix1 = LinearSpace.real.matrix(2, 3)(
             1.0, 0.0, 3.0,
             4.0, 6.0, 2.0
         )
         val matrix2 = 20.0 - (10.0 + (5.0 * matrix1))
         //val matrix2 = 10.0 + (5.0 * matrix1)
-        val expectedResult = Matrix.build(2, 3)(
+        val expectedResult = LinearSpace.real.matrix(2, 3)(
             5.0, 10.0, -5.0,
             -10.0, -20.0, 0.0
         )
@@ -76,15 +76,15 @@ internal class RealMatrixTest {
 
     @Test
     fun testSquareAndPower() {
-        val matrix1 = Matrix.build(2, 3)(
+        val matrix1 = LinearSpace.real.matrix(2, 3)(
             -1.0, 0.0, 3.0,
             4.0, -6.0, -2.0
         )
-        val matrix2 = Matrix.build(2, 3)(
+        val matrix2 = LinearSpace.real.matrix(2, 3)(
             1.0, 0.0, 9.0,
             16.0, 36.0, 4.0
         )
-        val matrix3 = Matrix.build(2, 3)(
+        val matrix3 = LinearSpace.real.matrix(2, 3)(
             -1.0, 0.0, 27.0,
             64.0, -216.0, -8.0
         )
@@ -95,16 +95,16 @@ internal class RealMatrixTest {
     @OptIn(UnstableKMathAPI::class)
     @Test
     fun testTwoMatrixOperations() {
-        val matrix1 = Matrix.build(2, 3)(
+        val matrix1 = LinearSpace.real.matrix(2, 3)(
             -1.0, 0.0, 3.0,
             4.0, -6.0, 7.0
         )
-        val matrix2 = Matrix.build(2, 3)(
+        val matrix2 = LinearSpace.real.matrix(2, 3)(
             1.0, 0.0, 3.0,
             4.0, 6.0, -2.0
         )
         val result = matrix1 * matrix2 + matrix1 - matrix2
-        val expectedResult = Matrix.build(2, 3)(
+        val expectedResult = LinearSpace.real.matrix(2, 3)(
             -3.0, 0.0, 9.0,
             16.0, -48.0, -5.0
         )
@@ -113,16 +113,16 @@ internal class RealMatrixTest {
 
     @Test
     fun testColumnOperations() {
-        val matrix1 = Matrix.build(2, 4)(
+        val matrix1 = LinearSpace.real.matrix(2, 4)(
             -1.0, 0.0, 3.0, 15.0,
             4.0, -6.0, 7.0, -11.0
         )
-        val matrix2 = Matrix.build(2, 5)(
+        val matrix2 = LinearSpace.real.matrix(2, 5)(
             -1.0, 0.0, 3.0, 15.0, -1.0,
             4.0, -6.0, 7.0, -11.0, 4.0
         )
-        val col1 = Matrix.build(2, 1)(0.0, -6.0)
+        val col1 = LinearSpace.real.matrix(2, 1)(0.0, -6.0)
-        val cols1to2 = Matrix.build(2, 2)(
+        val cols1to2 = LinearSpace.real.matrix(2, 2)(
             0.0, 3.0,
             -6.0, 7.0
         )
@@ -147,7 +147,7 @@ internal class RealMatrixTest {
 
     @Test
     fun testAllElementOperations() {
-        val matrix1 = Matrix.build(2, 4)(
+        val matrix1 = LinearSpace.real.matrix(2, 4)(
             -1.0, 0.0, 3.0, 15.0,
             4.0, -6.0, 7.0, -11.0
         )

kmath-for-real/src/commonTest/kotlin/kaceince/kmath/real/RealVectorTest.kt

@@ -2,9 +2,7 @@ package kaceince.kmath.real
 
 import space.kscience.kmath.linear.LinearSpace
 import space.kscience.kmath.linear.asMatrix
-import space.kscience.kmath.linear.real
 import space.kscience.kmath.linear.transpose
-import space.kscience.kmath.operations.invoke
 import space.kscience.kmath.real.plus
 import space.kscience.kmath.structures.Buffer
 import kotlin.test.Test
@@ -32,7 +30,7 @@ internal class RealVectorTest {
         val vector2 =  Buffer.real(5) { 5 - it.toDouble() }
         val matrix1 = vector1.asMatrix()
         val matrix2 = vector2.asMatrix().transpose()
-        val product = LinearSpace.real { matrix1 dot matrix2 }
+        val product = LinearSpace.real.run { matrix1 dot matrix2 }
         assertEquals(5.0, product[1, 0])
         assertEquals(6.0, product[2, 2])
     }